{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### **Description**\n", "\n", "This recipe shows how to calculate the upper Ocean Heat Content (OHC) using monthly-mean outputs from the National Oceanography Centre Near-Present-Day global eORCA1 configuration of NEMO forced using JRA55-do from 1976-2024.\n", "\n", "For more details on this model configuration and the available outputs, users can explore the Near-Present-Day documentation [here](https://noc-msm.github.io/NOC_Near_Present_Day/).\n", "\n", "---\n", "### **Background**\n", "Ocean Heat Content (OHC) is routinely used to characterise the energy stored by the ocean and can be defined using the conservative temperature field at time $t$ as follows:\n", "\n", "$$H(\\lambda, \\phi, t) = c_p \\int_{z=-h}^{\\eta} \\rho(\\lambda, \\phi, z, t) \\ \\Theta (\\lambda, \\phi, z, t) \\ dz$$\n", "\n", "where $c_p$ is the specific heat capacity of seawater, $-h$ is the lower depth, $\\eta$ is the sea surface, $\\rho$ is the in-situ density and $\\Theta$ is the conservative temperature.\n", "\n", "For a Boussinesq ocean model, such as NEMO, the above equation can be simplified by replacing the in-situ density with the reference density $\\rho_o$ and a constant ocean model heat capacity $c_p^o$ in place of $c_p$. \n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# -- Import required packages -- #\n", "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "\n", "from nemo_cookbook import NEMODataTree\n", "\n", "xr.set_options(display_style=\"text\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Using Dask**\n", "\n", "**Optional: Connect Client to Dask Local Cluster to run analysis in parallel.**\n", "\n", "Note that, although using Dask is not strictly necessary for this simple example using eORCA1, if we wanted to generalise this recipe to eORCA025 or eORCA12 outputs, using Dask would be essential to avoid unnecessary slow calculations using only a single process." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# -- Initialise Dask Local Cluster -- #\n", "import os\n", "\n", "import dask\n", "from dask.distributed import Client, LocalCluster\n", "\n", "# Update temporary directory for Dask workers:\n", "dask.config.set({'temporary_directory': f\"{os.getcwd()}/dask_tmp\",\n", " 'local_directory': f\"{os.getcwd()}/dask_tmp\"\n", " })\n", "\n", "# Create Local Cluster:\n", "cluster = LocalCluster(n_workers=4, threads_per_worker=3, memory_limit='5GB')\n", "client = Client(cluster)\n", "client" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Accessing NEMO Model Data**\n", "**Let's begin by loading the grid variables for our eORCA1 NEMO model from the [JASMIN Object Store](https://help.jasmin.ac.uk/docs/short-term-project-storage/using-the-jasmin-object-store/)**. \n", "\n", "**Alternatively, you can replace the ``domain_filepath`` below with a file path to your domain_cfg.nc file and read this with xarray's ``open_dataset()`` function.**" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

<xarray.Dataset> Size: 709MB\n",
       "Dimensions:       (y: 331, x: 360, nav_lev: 75)\n",
       "Coordinates:\n",
       "  * nav_lev       (nav_lev) int64 600B 0 1 2 3 4 5 6 7 ... 68 69 70 71 72 73 74\n",
       "  * x             (x) int64 3kB 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359\n",
       "  * y             (y) int64 3kB 0 1 2 3 4 5 6 7 ... 324 325 326 327 328 329 330\n",
       "Data variables: (12/49)\n",
       "    atlmsk        (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    bathy_metry   (y, x) float32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    bottom_level  (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    e1f           (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    e1t           (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    e1u           (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    ...            ...\n",
       "    top_level     (y, x) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    umask         (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "    umaskutil     (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    vmask         (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "    vmaskutil     (y, x) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    wmask         (nav_lev, y, x) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "Attributes:\n",
       "    CfgName:    UNKNOWN\n",
       "    CfgIndex:   -999\n",
       "    Iperio:     1\n",
       "    Jperio:     0\n",
       "    NFold:      1\n",
       "    NFtype:     F\n",
       "    VertCoord:  zps\n",
       "    IsfCav:     0\n",
       "    file_name:  mesh_mask.nc\n",
       "    TimeStamp:  01/03/2025 22:19:49 -0000

" ], "text/plain": [ " Size: 709MB\n", "Dimensions: (y: 331, x: 360, nav_lev: 75)\n", "Coordinates:\n", " * nav_lev (nav_lev) int64 600B 0 1 2 3 4 5 6 7 ... 68 69 70 71 72 73 74\n", " * x (x) int64 3kB 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359\n", " * y (y) int64 3kB 0 1 2 3 4 5 6 7 ... 324 325 326 327 328 329 330\n", "Data variables: (12/49)\n", " atlmsk (y, x) int8 119kB dask.array\n", " bathy_metry (y, x) float32 477kB dask.array\n", " bottom_level (y, x) int32 477kB dask.array\n", " e1f (y, x) float64 953kB dask.array\n", " e1t (y, x) float64 953kB dask.array\n", " e1u (y, x) float64 953kB dask.array\n", " ... ...\n", " top_level (y, x) int32 477kB dask.array\n", " umask (nav_lev, y, x) int8 9MB dask.array\n", " umaskutil (y, x) int8 119kB dask.array\n", " vmask (nav_lev, y, x) int8 9MB dask.array\n", " vmaskutil (y, x) int8 119kB dask.array\n", " wmask (nav_lev, y, x) int8 9MB dask.array\n", "Attributes:\n", " CfgName: UNKNOWN\n", " CfgIndex: -999\n", " Iperio: 1\n", " Jperio: 0\n", " NFold: 1\n", " NFtype: F\n", " VertCoord: zps\n", " IsfCav: 0\n", " file_name: mesh_mask.nc\n", " TimeStamp: 01/03/2025 22:19:49 -0000" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define directory path to ancillary files:\n", "domain_filepath = \"https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/domain_cfg\"\n", "\n", "# Open eORCA1 NEMO model domain_cfg:\n", "ds_domain = xr.open_zarr(domain_filepath, consolidated=True, chunks={})\n", "\n", "ds_domain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Next, we need to import the sea water conservative temperature and vertical grid scale factors stored at T-points in a single dataset.**" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

<xarray.Dataset> Size: 467GB\n",
       "Dimensions:                (time_counter: 577, y: 331, x: 360, deptht: 75,\n",
       "                            axis_nbounds: 2)\n",
       "Coordinates:\n",
       "  * deptht                 (deptht) float32 300B 0.5058 1.556 ... 5.902e+03\n",
       "    nav_lat                (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    nav_lon                (y, x) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "    time_centered          (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>\n",
       "  * time_counter           (time_counter) datetime64[ns] 5kB 1976-01-16T12:00...\n",
       "Dimensions without coordinates: y, x, axis_nbounds\n",
       "Data variables: (12/74)\n",
       "    berg_latent_heat_flux  (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "    deptht_bounds          (deptht, axis_nbounds) float32 600B dask.array<chunksize=(75, 2), meta=np.ndarray>\n",
       "    e3t                    (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 75, 331, 360), meta=np.ndarray>\n",
       "    empmr                  (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "    evs                    (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "    ficeberg               (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "    ...                     ...\n",
       "    ttrd_qns_li            (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "    vohfcisf               (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 75, 331, 360), meta=np.ndarray>\n",
       "    vohflisf               (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 75, 331, 360), meta=np.ndarray>\n",
       "    vowflisf               (time_counter, deptht, y, x) float32 21GB dask.array<chunksize=(1, 75, 331, 360), meta=np.ndarray>\n",
       "    zos                    (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "    zossq                  (time_counter, y, x) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "Attributes:\n",
       "    name:         OUTPUT/eORCA1_1m_grid_T\n",
       "    description:  ocean T grid variables\n",
       "    title:        ocean T grid variables\n",
       "    Conventions:  CF-1.6\n",
       "    timeStamp:    2024-Dec-09 10:18:50 GMT\n",
       "    uuid:         d4f9020a-a621-4a95-81da-0b477f67ed38

" ], "text/plain": [ " Size: 467GB\n", "Dimensions: (time_counter: 577, y: 331, x: 360, deptht: 75,\n", " axis_nbounds: 2)\n", "Coordinates:\n", " * deptht (deptht) float32 300B 0.5058 1.556 ... 5.902e+03\n", " nav_lat (y, x) float64 953kB dask.array\n", " nav_lon (y, x) float64 953kB dask.array\n", " time_centered (time_counter) datetime64[ns] 5kB dask.array\n", " * time_counter (time_counter) datetime64[ns] 5kB 1976-01-16T12:00...\n", "Dimensions without coordinates: y, x, axis_nbounds\n", "Data variables: (12/74)\n", " berg_latent_heat_flux (time_counter, y, x) float32 275MB dask.array\n", " deptht_bounds (deptht, axis_nbounds) float32 600B dask.array\n", " e3t (time_counter, deptht, y, x) float32 21GB dask.array\n", " empmr (time_counter, y, x) float32 275MB dask.array\n", " evs (time_counter, y, x) float32 275MB dask.array\n", " ficeberg (time_counter, y, x) float32 275MB dask.array\n", " ... ...\n", " ttrd_qns_li (time_counter, y, x) float32 275MB dask.array\n", " vohfcisf (time_counter, deptht, y, x) float32 21GB dask.array\n", " vohflisf (time_counter, deptht, y, x) float32 21GB dask.array\n", " vowflisf (time_counter, deptht, y, x) float32 21GB dask.array\n", " zos (time_counter, y, x) float32 275MB dask.array\n", " zossq (time_counter, y, x) float32 275MB dask.array\n", "Attributes:\n", " name: OUTPUT/eORCA1_1m_grid_T\n", " description: ocean T grid variables\n", " title: ocean T grid variables\n", " Conventions: CF-1.6\n", " timeStamp: 2024-Dec-09 10:18:50 GMT\n", " uuid: d4f9020a-a621-4a95-81da-0b477f67ed38" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define directory path to model output files:\n", "gridT_url = \"https://noc-msm-o.s3-ext.jc.rl.ac.uk/npd-eorca1-jra55v1/T1m\"\n", "\n", "# Construct NEMO model T-grid dataset, including conservative temperature (degC):\n", "# (Combine 75 vertical levels into single chunk for transforming vertical coordinates with dask)\n", "ds_gridT = xr.open_zarr(gridT_url, consolidated=True, chunks={'deptht': 75})\n", "\n", "ds_gridT" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Creating a NEMODataTree**\n", "\n", "**Next, let's create a NEMODataTree to store our domain and T-grid variables for the eORCA1 model.**" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

<xarray.DataTree 'NEMO model'>\n",
       "Group: /\n",
       "│   Dimensions:               (time_counter: 577, axis_nbounds: 2)\n",
       "│   Coordinates:\n",
       "│       time_centered         (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>\n",
       "│     * time_counter          (time_counter) datetime64[ns] 5kB 1976-01-16T12:00:...\n",
       "│   Dimensions without coordinates: axis_nbounds\n",
       "│   Data variables:\n",
       "│       time_centered_bounds  (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array<chunksize=(1, 2), meta=np.ndarray>\n",
       "│       time_counter_bounds   (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array<chunksize=(1, 2), meta=np.ndarray>\n",
       "│   Attributes:\n",
       "│       nftype:   F\n",
       "│       iperio:   True\n",
       "├── Group: /gridT\n",
       "│       Dimensions:                (time_counter: 577, j: 331, i: 360, k: 75,\n",
       "│                                   axis_nbounds: 2)\n",
       "│       Coordinates:\n",
       "│         * deptht                 (k) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03\n",
       "│           time_centered          (time_counter) datetime64[ns] 5kB dask.array<chunksize=(1,), meta=np.ndarray>\n",
       "│           gphit                  (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           glamt                  (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│         * k                      (k) int64 600B 1 2 3 4 5 6 7 ... 69 70 71 72 73 74 75\n",
       "│         * j                      (j) int64 3kB 1 2 3 4 5 6 ... 326 327 328 329 330 331\n",
       "│         * i                      (i) int64 3kB 1 2 3 4 5 6 ... 355 356 357 358 359 360\n",
       "│       Dimensions without coordinates: axis_nbounds\n",
       "│       Data variables: (12/80)\n",
       "│           berg_latent_heat_flux  (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "│           deptht_bounds          (k, axis_nbounds) float32 600B dask.array<chunksize=(75, 2), meta=np.ndarray>\n",
       "│           e3t                    (time_counter, k, j, i) float32 21GB dask.array<chunksize=(1, 75, 331, 360), meta=np.ndarray>\n",
       "│           empmr                  (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "│           evs                    (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "│           ficeberg               (time_counter, j, i) float32 275MB dask.array<chunksize=(1, 331, 360), meta=np.ndarray>\n",
       "│           ...                     ...\n",
       "│           e1t                    (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           e2t                    (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           top_level              (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           bottom_level           (j, i) int32 477kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           tmask                  (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "│           tmaskutil              (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│       Attributes:\n",
       "│           nftype:   F\n",
       "│           iperio:   True\n",
       "├── Group: /gridU\n",
       "│       Dimensions:       (j: 331, i: 360, k: 75)\n",
       "│       Coordinates:\n",
       "│           gphiu         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           glamu         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│         * k             (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75\n",
       "│         * j             (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331\n",
       "│         * i             (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5\n",
       "│       Data variables:\n",
       "│           e1u           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           e2u           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           umask         (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "│           umaskutil     (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│       Attributes:\n",
       "│           nftype:   F\n",
       "│           iperio:   True\n",
       "├── Group: /gridV\n",
       "│       Dimensions:       (j: 331, i: 360, k: 75)\n",
       "│       Coordinates:\n",
       "│           gphiv         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           glamv         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│         * k             (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75\n",
       "│         * j             (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5\n",
       "│         * i             (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360\n",
       "│       Data variables:\n",
       "│           e1v           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           e2v           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           vmask         (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "│           vmaskutil     (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│       Attributes:\n",
       "│           nftype:   F\n",
       "│           iperio:   True\n",
       "├── Group: /gridW\n",
       "│       Dimensions:       (j: 331, i: 360, k: 75)\n",
       "│       Coordinates:\n",
       "│           gphit         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           glamt         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│         * k             (k) float64 600B 0.5 1.5 2.5 3.5 4.5 ... 71.5 72.5 73.5 74.5\n",
       "│         * j             (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331\n",
       "│         * i             (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360\n",
       "│       Data variables:\n",
       "│           e1t           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           e2t           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│           wmask         (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "│           wmaskutil     (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "│       Attributes:\n",
       "│           nftype:   F\n",
       "│           iperio:   True\n",
       "└── Group: /gridF\n",
       "        Dimensions:       (j: 331, i: 360, k: 75)\n",
       "        Coordinates:\n",
       "            gphif         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "            glamf         (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "          * k             (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75\n",
       "          * j             (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5\n",
       "          * i             (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5\n",
       "        Data variables:\n",
       "            e1f           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "            e2f           (j, i) float64 953kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "            fmask         (k, j, i) int8 9MB dask.array<chunksize=(75, 331, 360), meta=np.ndarray>\n",
       "            fmaskutil     (j, i) int8 119kB dask.array<chunksize=(331, 360), meta=np.ndarray>\n",
       "        Attributes:\n",
       "            nftype:   F\n",
       "            iperio:   True

" ], "text/plain": [ "\n", "Group: /\n", "│ Dimensions: (time_counter: 577, axis_nbounds: 2)\n", "│ Coordinates:\n", "│ time_centered (time_counter) datetime64[ns] 5kB dask.array\n", "│ * time_counter (time_counter) datetime64[ns] 5kB 1976-01-16T12:00:...\n", "│ Dimensions without coordinates: axis_nbounds\n", "│ Data variables:\n", "│ time_centered_bounds (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array\n", "│ time_counter_bounds (time_counter, axis_nbounds) datetime64[ns] 9kB dask.array\n", "│ Attributes:\n", "│ nftype: F\n", "│ iperio: True\n", "├── Group: /gridT\n", "│ Dimensions: (time_counter: 577, j: 331, i: 360, k: 75,\n", "│ axis_nbounds: 2)\n", "│ Coordinates:\n", "│ * deptht (k) float32 300B 0.5058 1.556 ... 5.698e+03 5.902e+03\n", "│ time_centered (time_counter) datetime64[ns] 5kB dask.array\n", "│ gphit (j, i) float64 953kB dask.array\n", "│ glamt (j, i) float64 953kB dask.array\n", "│ * k (k) int64 600B 1 2 3 4 5 6 7 ... 69 70 71 72 73 74 75\n", "│ * j (j) int64 3kB 1 2 3 4 5 6 ... 326 327 328 329 330 331\n", "│ * i (i) int64 3kB 1 2 3 4 5 6 ... 355 356 357 358 359 360\n", "│ Dimensions without coordinates: axis_nbounds\n", "│ Data variables: (12/80)\n", "│ berg_latent_heat_flux (time_counter, j, i) float32 275MB dask.array\n", "│ deptht_bounds (k, axis_nbounds) float32 600B dask.array\n", "│ e3t (time_counter, k, j, i) float32 21GB dask.array\n", "│ empmr (time_counter, j, i) float32 275MB dask.array\n", "│ evs (time_counter, j, i) float32 275MB dask.array\n", "│ ficeberg (time_counter, j, i) float32 275MB dask.array\n", "│ ... ...\n", "│ e1t (j, i) float64 953kB dask.array\n", "│ e2t (j, i) float64 953kB dask.array\n", "│ top_level (j, i) int32 477kB dask.array\n", "│ bottom_level (j, i) int32 477kB dask.array\n", "│ tmask (k, j, i) int8 9MB dask.array\n", "│ tmaskutil (j, i) int8 119kB dask.array\n", "│ Attributes:\n", "│ nftype: F\n", "│ iperio: True\n", "├── Group: /gridU\n", "│ Dimensions: (j: 331, i: 360, k: 75)\n", "│ Coordinates:\n", "│ gphiu (j, i) float64 953kB dask.array\n", "│ glamu (j, i) float64 953kB dask.array\n", "│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75\n", "│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331\n", "│ * i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5\n", "│ Data variables:\n", "│ e1u (j, i) float64 953kB dask.array\n", "│ e2u (j, i) float64 953kB dask.array\n", "│ umask (k, j, i) int8 9MB dask.array\n", "│ umaskutil (j, i) int8 119kB dask.array\n", "│ Attributes:\n", "│ nftype: F\n", "│ iperio: True\n", "├── Group: /gridV\n", "│ Dimensions: (j: 331, i: 360, k: 75)\n", "│ Coordinates:\n", "│ gphiv (j, i) float64 953kB dask.array\n", "│ glamv (j, i) float64 953kB dask.array\n", "│ * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75\n", "│ * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5\n", "│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360\n", "│ Data variables:\n", "│ e1v (j, i) float64 953kB dask.array\n", "│ e2v (j, i) float64 953kB dask.array\n", "│ vmask (k, j, i) int8 9MB dask.array\n", "│ vmaskutil (j, i) int8 119kB dask.array\n", "│ Attributes:\n", "│ nftype: F\n", "│ iperio: True\n", "├── Group: /gridW\n", "│ Dimensions: (j: 331, i: 360, k: 75)\n", "│ Coordinates:\n", "│ gphit (j, i) float64 953kB dask.array\n", "│ glamt (j, i) float64 953kB dask.array\n", "│ * k (k) float64 600B 0.5 1.5 2.5 3.5 4.5 ... 71.5 72.5 73.5 74.5\n", "│ * j (j) int64 3kB 1 2 3 4 5 6 7 8 ... 325 326 327 328 329 330 331\n", "│ * i (i) int64 3kB 1 2 3 4 5 6 7 8 ... 354 355 356 357 358 359 360\n", "│ Data variables:\n", "│ e1t (j, i) float64 953kB dask.array\n", "│ e2t (j, i) float64 953kB dask.array\n", "│ wmask (k, j, i) int8 9MB dask.array\n", "│ wmaskutil (j, i) int8 119kB dask.array\n", "│ Attributes:\n", "│ nftype: F\n", "│ iperio: True\n", "└── Group: /gridF\n", " Dimensions: (j: 331, i: 360, k: 75)\n", " Coordinates:\n", " gphif (j, i) float64 953kB dask.array\n", " glamf (j, i) float64 953kB dask.array\n", " * k (k) int64 600B 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75\n", " * j (j) float64 3kB 1.5 2.5 3.5 4.5 ... 328.5 329.5 330.5 331.5\n", " * i (i) float64 3kB 1.5 2.5 3.5 4.5 ... 357.5 358.5 359.5 360.5\n", " Data variables:\n", " e1f (j, i) float64 953kB dask.array\n", " e2f (j, i) float64 953kB dask.array\n", " fmask (k, j, i) int8 9MB dask.array\n", " fmaskutil (j, i) int8 119kB dask.array\n", " Attributes:\n", " nftype: F\n", " iperio: True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define dictionary of grid datasets defining eORCA1 parent model domain with no child/grand-child nests:\n", "datasets = {\"parent\": {\"domain\": ds_domain, \"gridT\": ds_gridT}}\n", "\n", "# Initialise a new NEMODataTree whose parent domain is zonally periodic & north-folding on F-points:\n", "nemo = NEMODataTree.from_datasets(datasets=datasets, iperio=True, nftype=\"F\", read_mask=True)\n", "\n", "nemo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Calculating Upper Ocean Heat Content**\n", "\n", "**Now we have constructed our `NEMODataTree`, let's calculate the Ocean Heat Content of the upper 500 m of the Pacific Ocean.**\n", "\n", "$$H_{upper}(t) = \\rho_o c_p^o \\sum_{z^{*}=0}^{z=500} \\sum_{j} \\sum_{i} (e_{1t}(j, i) \\ . \\ e_{2t}(j, i) \\ . \\ e_{3t}(t, z, j, i) \\ . \\ \\Theta(t, z, j, i))$$\n", "\n", "In this example, our eORCA1 model uses $z^{*}$ vertical coordinates with time-evolving vertical grid scale factors, hence $e_{3t}(t, z, j, i)$ is a function of time and space.\n", "\n", "Let's start by defining our Pacific Ocean basin mask:\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

<xarray.DataArray 'pacmsk' (j: 331, i: 360)> Size: 119kB\n",
       "dask.array<astype, shape=(331, 360), dtype=bool, chunksize=(331, 360), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * i        (i) int64 3kB 1 2 3 4 5 6 7 8 9 ... 353 354 355 356 357 358 359 360\n",
       "  * j        (j) int64 3kB 1 2 3 4 5 6 7 8 9 ... 324 325 326 327 328 329 330 331

" ], "text/plain": [ " Size: 119kB\n", "dask.array\n", "Coordinates:\n", " * i (i) int64 3kB 1 2 3 4 5 6 7 8 9 ... 353 354 355 356 357 358 359 360\n", " * j (j) int64 3kB 1 2 3 4 5 6 7 8 9 ... 324 325 326 327 328 329 330 331" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define Pacific Ocean basin mask:\n", "pacmask = ds_domain['pacmsk'].rename({\"x\":\"i\", \"y\":\"j\"}).astype(bool)\n", "# Assign (i,j) coordinates of T-points:\n", "pacmask['i'] = pacmask['i'] + 1\n", "pacmask['j'] = pacmask['j'] + 1\n", "\n", "pacmask" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of estimating the upper 500 m ocean heat content by integrating over the `Nk` uppermost vertical levels, such that `deptht[Nk]` is closest to 500 m depth, we can use the `.depth_integral()` method to perform vertical integration between two depth surfaces.\n", "\n", "Here, we:\n", "\n", "- Select the time period between 2020-2024.\n", "- Apply the Pacific Ocean mask. \n", "- Integrate first with respect to depth (between 0 - 500 m).\n", "- Integrate over the NEMO model T-grid ($i$, $j$) dimensions." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/z2/j_dr250s42x34hk63_rp4bm80000gq/T/ipykernel_36206/1206376866.py:7: RuntimeWarning: Indexing with a boolean dask array is not allowed. Mask will be computed first using .compute(). This may result in high memory usage for large masks.\n", " .apply_mask(pacmask, drop=True)\n" ] }, { "data": { "text/html": [ "\n", "

\n", " \n", "

\n", " \n", " NEMODataTree 'NEMO model'\n", "

\n", "\n", "

\n", " NEMODataArray 'integral_ij(integral_z(thetao_con))'\n", "

├─ Domain: '.'

\n", "

├─ Grid Path: 'gridT'

\n", "

└─ Grid Type: 'T'

\n", "

\n", " \n", "

<xarray.DataArray 'integral_ij(integral_z(thetao_con))' (time_counter: 49)> Size: 392B\n",
       "dask.array<mul, shape=(49,), dtype=float64, chunksize=(1,), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * time_counter   (time_counter) datetime64[ns] 392B 2020-01-16T12:00:00 ......\n",
       "    time_centered  (time_counter) datetime64[ns] 392B 2020-01-16T12:00:00 ......

\n", " \n", " " ], "text/plain": [ "\n", " \n", "\n", " Size: 392B\n", "dask.array\n", "Coordinates:\n", " * time_counter (time_counter) datetime64[ns] 392B 2020-01-16T12:00:00 ......\n", " time_centered (time_counter) datetime64[ns] 392B 2020-01-16T12:00:00 ......" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Calcuate upper 500m ocean heat content [Zeta Joules = 10^21 Joules] for the entire Pacific Ocean:\n", "rho0 = 1025 # Reference density [kg m-3]\n", "cp0 = 3991.86 # Reference specific heat capacity [J kg-1 C-1]\n", "\n", "H_upper = rho0 * cp0 * (nemo[\"gridT/thetao_con\"]\n", " .sel(time_counter=slice(\"2020-01\", \"2024-01\"))\n", " .apply_mask(pacmask, drop=True)\n", " .depth_integral(limits=(0, 500))\n", " .integral(dims=[\"i\", \"j\"])\n", " )\n", "\n", "H_upper" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Notice that the output above contains dask arrays, so we haven't actually computed the Upper Ocean Heat Content yet. To do this, we need to call the ``.compute()`` method:**" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "H_upper = H_upper.compute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Visualising the evolution of Pacific Upper Ocean Heat Content**\n", "\n", "**Finally, let's visualise the results by plotting the time-series of Pacific Upper Ocean Heat Content:**" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Upper Ocean Heat Content (ZJ)')" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAGaCAYAAADuEJ7dAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAoIdJREFUeJztnQm4VWPbxx9NhpTmOYrMYyjzUCkNypgpiZBMUa8pXhRvhQ8v8ZppUCTSgERSkcosSRlCpFTSJFFyvuv3fO+zvnX22fucvc/Zw1pr/3/XtWpPZ++1117D/7mf/33f2xQUFBQYIYQQQggh8ohyuV4BIYQQQgghso1EsBBCCCGEyDskgoUQQgghRN4hESyEEEIIIfIOiWAhhBBCCJF3SAQLIYQQQoi8QyJYCCGEEELkHRLBQgghhBAi75AIFkIIIYQQeYdEsBB5xB9//GEefvhh06pVK1OrVi1TsWJF+//xxx9v/vOf/9jnY2nSpInZZpttCi3bbbedadasmbn88svNsmXL4n4WzShfeuklc+qpp5pGjRqZbbfd1tSsWdMceOCB5h//+If54osvivzN77//bqpWrep9Trly5cz3338f9/2nTp1qzjvvPLPHHnvY17m/mTFjRtLb44ILLvD+bsCAAd7j3I79zhUqVDB16tQxHTt2NG+88Uax77v77rsX+tuZM2eW+PnxltjvXtxr+Q1LyzfffGMuu+wyc+ihh5pKlSrF3SZlId735DerXr263Reff/55k2nc57I/+1m5cqU5//zzTYMGDexvzGtOOeUUu+3TsW0T8c4775hzzz3Xrg/HU5UqVcy+++5rrr76avPtt9+m/fOEEHGgbbIQIvr8+OOPBfvvvz9t0hMuPL906dJCf7fLLrsU+zdNmjQp+O233wr9zdq1awvatGlT7N9dffXVRdZx5MiRRV43YMCAuN+Hv4/3vtOnT096m/To0cP7u9tuu817nNvFrfs222xT8PLLL8d9z7fffrvI6y+44IISPz/e8t133xV6fXGvPe644wpKy/jx4+O+p3+blIWSvifLjTfeWJBJ3OewP/s5++yzi6zLySefbLd9OrZtPPr161fstth+++0LxowZk9bPFEIURZFgIfIANMBpp51m5s+fb++3bNnSfPLJJ+bPP/+0/7do0cI+zvO87v80Q1GmT59u/vrrLzN79mwbuQIiZhMnTiz0ujPPPNNMmzbN3q5Xr54ZM2aMWbt2rdm0aZP56KOPTL9+/cz2229f5P2HDx/u3SYCByNGjIi7Pocccoi54447zOuvv26jy5miR48e9vN//vln0759e/sY94cOHRr39fG+w4svvmg2btyY8DN22WUX+56xS2zU0nHbbbcVeW0qEfBYGjZsaK6//nrzwgsvmIsuushkErfuGzZsKBRpvvvuu80PP/yQsc912yk2us7+6Pjuu+/sayZMmGC3fTq2bSwPPvigue++++ztHXbYwTzzzDPmt99+MytWrDA33HCDfZzjhOj0xx9/nLbPFULEIY4wFkKEhG+//bagV69eBU2bNi2oVKlSQZUqVQqOOeaYgrFjxxZ63YQJE7woE6/76aefikSJK1as6L1m4sSJcSPB/ijrqaee6j0+ePBg7/E33nijUMR0zpw5cdd9y5Ythe4vWbKkoFy5cvbv+D4nnnii9z4zZswodjscdthhGYsE8xoH28U9vsceexR5v40bN9rfgOcrV65c0L17d+/1w4cPT/j5sdHJRCQbofX/Nu+++26h5y6++GLvuXjR7BtuuKHEzyEymsr2TrSd//77b297sbzwwgtelL9FixYFderUsfsrkVG2N4+vWrWqyPvPnTu34Kyzzipo2LCh3Y+rV69e0LJlS+/9/NvObethw4YljMSyjsVFgn/99deCm266yc6c7LDDDgXbbbddwW677VZw6aWXFrsd2Odr167tve8999xT5DXt27f3nud3jD1Or7rqKrst+Ez2sb333rvI7/Tpp58WnHvuuYW2B8fTm2++Weh1y5Yts6/bZ599CmrUqFFQoUIF+3sccsghBffee2+hYzR2e0yZMqXg8MMPt+vBNu3fv3/B5s2bi/3+QgQNRYKFCCkffPCBjYA+/vjjNoK1efNmG13Da0gktn///t5rX3vtNe92hw4drP/RD55dF+WEKVOmlPj5/ugsXlmHPyqM3/Pwww+P+/f4L/0Q8f3777/t7bPOOst+h3jR1VyS6Ds7xo0bZ38D6Ny5s/XCJvMd8FU7jzZRWbzOX375ZcLXE4Umks6y//77m0GDBtmovqNXr17ebSKNDl5DtNf95uwLQYRtxf6NX5f9msjoV199ZR544AHrz92yZYv32ieeeMIceeSR1lf8008/2efWrFlj3n//fTNr1qy0rxuRZI67wYMH25kTfOx46RcvXmxnPIrjww8/NKtWrfLu9+zZs8hrLrnkEu823nN3TBCx5rcmksy24DOZXVi4cKGNXDsmTZpkZ3aeffbZQtuDGZO2bduaRx991Hst0Wdehz//119/tbM87L98Fr79K664Iu73YPYIb/zcuXPteixZssQMGTLE3HPPPSluTSFyi0SwECGFCygXrGrVqpk333zTXoyYTj7mmGPs83fddZf5/PPP7W3/NPOuu+4a9/1222037zYXtURs3brVXvz4TMAWcfLJJ3vPI8gdJPoky8iRI73biGAS6hCFydgJsgGC4ZFHHvHud+/evchr/EKX73DccceZunXr2vskxyVK8kOorF692ooQBPHo0aNtklqi6XBEDb83C7/xP//5Tytw+G2gXbt2npUCcYiQdAJp3bp19ja2h/Lly5tcwe85cOBAb9BAohw2HWA7L1q0yFpo2DZLly71BmkLFizwBmlsqz59+nhC8aabbjLLly+3f4eAPOKIIxJ+PgMUBjVYURzO/lBcQiCf9+OPP9rbDPD4jfguCEmEY3H4j0OSAlmKOw55X/YLuPDCC+33gk6dOtntw/N8PtYJYLBw8cUX223G789AgoEPA6o999zTfjesSL/88ot9PYNhBm4c77wXr/3ss8/sAAmeeuop7zP9rF+/3vTt29cKZ78AZyArRKjIdShaCJE6X3/9dYmJRv7pVv8Ua9++fUtMNOvYsWNSiXH77rtvwezZswu9T4cOHbzn+/Tpk9T38SeT+W0GrEdxdoJs2CFilx133LFg0KBBRd7Lb+eoWrVqwaZNm+zjl19+ecIkvxEjRljryg8//GBf/9FHHxUcfPDB3uvbtm1b6PUkj7G9mY4n+fCZZ56xdgH3+tGjR3uvveOOO7zHx40bZx876aST7P3y5cvbz4xHMnaITCTGXX/99YUS9U444YSCWrVq2XWNfe2dd95pX/fkk096jx1//PHFrkOsHSLe/u0nnh2C3wjLgHv8+++/T2k7kOzm/haLQjw++eSTQt915cqVBd988413H7vChg0b4v7t1KlTkzovvPjii/b12B04R2A9YZ/FvhT7WqwmsdsDS4ffKlGzZk37+LbbbpvS9hAi1ygSLEQIISqZDC7i4492JSq/5H/c//riIHrkn5qOjWQRtUsGfwSV6N2nn35ql4MPPjjua3IJ0VYSmWLx2zkOO+wwG6njO+y3336FXuO3VBDB69q1q2ncuLEtk8X39SfckYDohylntg8RxJ122snaJlwUMPb1zBS4SC+WCKbhXQQVGwSfmUtIGuQ7HHvssWbUqFF25sJF/ZkFYKaB/ddFt/0Q8QSSFR1YBTKNi9a7GZBkjxOH//VE8+NFWf3HYeXKlW1ZQf/3JMK74447puW8QDT32muvtRFjorvxElDdto4tAei3M7Ge4LfkCBEGJIKFCCFuih322muvuJUFWPAtgt/7iRBiytgP3kE8gw6/Pzi2OgTT1//617/sfab3u3TpYv/e4bdGvPXWW9abGQ8nJvBUOp+qE4rNmze3i/uckuwEmYTqEKwr/lK2O6IAMfrQQw8ltHNQw9h9B2op+60ib7/9tr3tBHNxYBFwxBODxb2eqW6mzWHy5Mm2DrTb5n7PcLZx1SH4/ohAftdu3bp5z2MFcVx33XWeOGMaPxYqjzic9SeTIEid+OM4SLWaBRaX2rVre/effvrpIq/BguDA1sJv6v+eHAOJrEH+88KJJ54Y95zAdr/00kvtaxh8OKjpjYjlNf7BZzycTSm2CooQYUMiWIgQQqMKF2Ek4kg0B2FLVJZIEg0xDjjgAM/bS5IWF2DgQnf66aebefPm2dfzP/edb5SkGl6fCKJQN998syd28ZhSXstxwgkn2AswcEElqofIRczgYSWpBkGDjzU2maw4eC+/5xDxTESLxYk7tz7u8XRBRPWoo44qlFTE+rvPIBmRhhPJ4CLaCCg8sERp8bbyu+DvpFmCA0+xg0gpPl5EI5Foticixp/45n+9X+zy25I8B/g9SWryw37gtpk/8uffxv6IP8lppWlOkgz+CCMlxBBcbN94flMGd0TQ3QDt1ltvtdFQtg33092Eg8/yb7tzzjnHRvtd4p5/0Jbou91yyy2FBgQkprGdidRzXDFgAb63O0aYXeF4Bo4VmmzweXwuHl5Xco191IlsPNEkqhG9Zt/iPME+xLnDvz4OItscRyTacowKkRfk2o8hhCgd77//vvXxJdtsIdlmGbzOT6ISaV999ZVXVg0v4XvvvVeqZhmtW7f2HnvggQeKfM/XXnvNe57SaZTVSqahRTKnN38Js4EDB5ZYIg3wnrrnrrjiCvtYz549i/VcL1q0qJCnmOYifo9lvAXfJX8Xb53iLTR4cNvGsXXr1oLGjRsXet2tt95aZP34XUvalv7fPl0l0kryzfoXvOLx3uPxxx/3vNiJ9rF0eYIBH3DsNnXLTjvtVJAM7CPFbWvKjj333HOF/ubDDz8sqFatWtzXH3jggd7rJk2aVMgnXtxx0bt37yLPUfKtUaNGRX7f4krGJdqGQgQdRYKFCClEbIkCMd1OdIe2xERp8evhMyXi6C+FRgQQawLT+EQMa9SoYSNB/I8vk9JLPO8yw0uCz3FT/WgMf2Y8Xk8iUUSAiRizHrTjxctKRAsvIlFKsuxdJJHn/dPi/ilh51/12wnSAeWl/FPdyXDvvfd607+PPfaYrQrgt3PEK3tFZj5lvIAoLr5Xprgp+UVkEZ8n5c6INPLaa665xv623PZXmyDizu/OtDe/HduZaiCUCWM6O3Zamql0f/ML7lM9IMjwPYm40w6bfZptwPcj6hoPSorhhebvKC9HBJWKKUTZjz766LSvH75eor9UomA2xv1uRGtZh2QgcktE/+yzz7b7Nvs+Ue+9997bXHXVVdbawXOxzWHYJ3jebRv+BjuU34LELA4lzvCK77zzznZ7sJ/w3jzmj46zL7OvcXzyHfCbY+Xx+/qFiDLboIRzvRJCCJEtmPJlypnayX57AwIjG8lVQgghgkHhavVCCBFxiMT6I2cumigBLIQQ+YVEsBAi78A2QCIQwpemCX7LgBBCiPxAdgghhBBCCJF3KDFOCCGEEELkHRLBQgSYAQMGePVYmbb3Q9F891wUi9Xzff3fL3aJbZxBju+TTz5pqwJQJYOFzm00H4iX/7ty5UrTp08fW5mBTPs6derYqhokyOWCl19+2dZUpsIB60NGP00LqB3r7xgW5d+Z/T0Wfp9M1SROxP3332/XJd76lASNLKj60apVK1uzl8oPVF+ghu+dd96Z1vrVxTFhwgTvO2S7yQy/k/tsKmkIEVTkCRZCRILevXvbQv9+KPnGQhMKOqY5aCxCyTK/OKBZAaXLXn31VVs5IrbxRKagCQXtj8eOHVvocRpc0LSAhSoWNBWhSYXIjgh2jWZSEcLz58+3SZeU8vPD/sZCKTdKmp1yyikmGyLYNRhhv2EwkU0RPHDgQHubzz3ooIOy9tlCpIIiwUKInECXrGTrssZr/+q/qNPy2Qlg6srSvY3F1dmlg96bb77pvZ46xU4AUyeVrlp07qKOLl24XKvkbHDjjTd6AphoNBFh1oGOcqwH/PrrrzZKvHTp0qysk0gd9iE6JToBzIwENa35LeliSP3d0047LZKzNkKEllx36xBCJKa4zmWxHcfiPU5np9dff72gZcuWtgtV3bp1bRet33//PW63MD5j9OjRtnPctttuaztj0UmNzmN+NmzYUDBgwAD7OjpM8d777bdfwZAhQwr+/PPPQq/1d+qaPXt2wbHHHmv/JrbrVKIuY7EdvuJxxhlneJ8zfPhw73Fuu8fPPPNM+9iaNWsKKlSo4HXm2rRpU9xucJMnT467PenIdcABB9jtw/+vvPKK3T7/8z//Y9eVjnBHH3207fBVEitXrizU3evFF18s9Dwd4PbZZ5+E3ei++OIL262uSZMm9n3oIEj3sP/85z+FXjdz5syCU045xf7+dPmjG91pp51WZB0///xzuy3pzkZ3svLly9suaEcddVTBU089VagjXex+M2rUKLs92KbNmjUruO+++4p0sCtNN7lEHQsz8b2GDRuWVKe1eNx0003e61ifdevWxX3dli1bvNscK+w3Bx98cEHlypXtPsW2u+qqqwqWL19e6O/8Xfo4js4///yCGjVq2O/Rvn37gm+++ca+rqROhP5tOGHChIJ27drZ9+GYaNCgge2iSDfIRL/RxIkT7frVq1cv7r5e3GezfYUIEhLBQkRYBHNx44IfezHq2LFjXDGDiIh38brsssu8169evbqQMItdELl+Iexvx7r99tsnbL0ai7vwIm5q1qzpXaS7detWqJ0w+NvYfvTRR97j3HaP77zzzvaxadOmeY8h4v1wcY9tLxy7PWNb9LJeCK/Y7VCrVq2EQshBa1z3+urVqxcZbMC9997rvWbvvff2Hn/11VetaErUQtnx8MMP27bW8V7Htn355Ze9144fP75YETN48OC4+w3rHu/1zz77bLHfvywiOBPfqywieN999/Ved/vtt5f4vRl8ISATfRZC2gnbWBEcb3uzb/z1119Ji+Abbrgh4WsQtx988EHc3yjeZ/v39eI+WyJYBA2JYCEiLIJZiORygSJ6xMXKPT5lypQiYoaFyBiRXr/IQmwsXLjQvv7KK6/0XvvQQw8VrF+/vmDt2rUFffr0KfS4w//erVq1slGmjRs3FixYsKDY7+6/8Ma7SPvFrj+aunjxYu9xRIR7nO8CCDP32DHHHFPoM//5z396z1166aVxtyfRbr6zXzCzPPDAA3Y7d+nSxXsMkVscd911l/fa5s2bx30N0Tr/QMIJqDp16niPEw3+/vvv7e82a9asgscee8y+bunSpd5vSLSR35ABCpE7N+CpX7++F51ke/G7//TTTwV//PGH/Rz2Gz6X1zIYcVHT2P2GyC/f/8EHH/QeO/HEE4v9/iX9zokEXCa/V6zwThb3Xizjxo0r8fV33313od+e4+KXX34pOO+88+IOVv0imFkXBoJsB8Sve3zOnDlxt2tsBB2B654jisy+w/ZjgOiOpUMPPTTuezVs2LDg/ffft+vqnznx7+v+85aErwgy8gQLEWHISr/llltM1apVzRFHHGE7ozneeOONIq/nNT179rSVFTp27Gg9jICWxdMI48eP915/5ZVX2veuVq2aGTp0qPf4lClT4q4PiTq777672WGHHcw+++xT7Lq3bt3aemXxxuKr/Oijj2y1BNf1DS9tPPyVIPy343kxY6tGlPT6+vXrm+uvv9422ujQoYP3OBUdrrrqKrst2G6OkrLyk+laH+817777rq1uAbvttpv1Q+Od5nejCkGvXr3scyT4/fnnn/Y2yYF77723rTxx6KGH2kRAIGFr3rx59na9evXM3Llz7XerWbOm/Z1IIHT+bXyv7nP9NG/e3Pqs+f7+KiaZqkqQre9VWpLx/U6cONG7TfIdxwXrxnHk/p5j9I8//ijyt4MGDbJ+d/a7VPY3f9Kc/1h1FVLatGljEzLhww8/jFvJgv2/RYsWdl1PP/30lD9biCAhESxEgNluu+0SJpJRiine6/zsvPPOhS7ICCVHvIu+//lEr1+xYkWJ6x3v4km5qMaNG5tkOf/8823JMv6G74cA9gttMu0ddevW9W6vXbvWu01CUuxrEr020ev97LrrrjZ5DrbffvsiZbwAMeGIJ2ASbV9ExN9//13kNd9++22R1/tLpjGYKF++fNz3T+a38v9eZ599trnjjjtsmTj2r3gCnAFJLIhQR+XKlZP+/rHcdtttRRIgY/fJbH6vVGDfcCxYsKDE1/u/g/87Vq9e3Q4mgORMkiLTvb2T3X4MDtL92UIECYlgIQKM/+IYW7/WRbliX+fnxx9/LHTBd2WfXCWCWPzPJ3q9E4eIvmXLlsWt3OAXqA6ib8kSTwzG4sQoHH744d7tzz//vFDJKgc1g+GQQw4xFStWtLepIOG/eMd7vZ8KFeJXlUz0eEkQeaOOLKxZs6ZQhA7Ylk8//bR3v3379l5k07Fw4cKE28sv5C+99NK4vxV/S1UDBgSvvPKKJ+RnzZply7fxmho1ahT7Pdz2hGxUP8j09yrNd6A0muOhhx4yGzZsiPs6V3XE/x38xxn7wfr16739Kt46JrO9i/sO/s8eMmRIwu3nqquk87OFCBISwUIEmLZt23oRxy+//NLWwqXU16hRowrZAbp06RL373/66Sc7dcpFlengJ554wnuuXbt2RV4/Z84cM3z4cGs3YMr5pZde8i5qrAtQqgu4UFLCCxGGqCA6SZ1dhNozzzxTpu+NBYISU7wPQpupb6a9r776au81/jq+F110UaGLOuJ28eLF5q677vIed1YQmlC4aVwEcP/+/a3wGDNmjHnnnXe8QYX7vpmE6DiWEscVV1xhJk+ebL8vvx3fy0UVWe9+/frZ21ge3KCE74oQZJsR5Xzvvfe8cnFM/7vI9LBhw8zIkSNttJuoJ00M/vnPf1pbgBNcTrwwwMDyweuIzsaLRuaSTH8vpvodyTZ74LfBLuMirQhwbCv8lhx/HLfsd9Shjj1mb7/9dvs7sj6U7HMDV47RRLM8JeH/Dgyg/QMlf53iu+++2w4S2Hc47jlPcJw5K1RZP5tBabbKDQqRMrk2JQshisefaBRvoezTr7/+6r3en8hFkhCZ8rF/06FDh7gJTlRfiPcZvXv3LlQdwp8JH2/xJ8O4x5IpdRbvO8Rb+F6xFSJ69eqV8PWXX355odcuW7bMlhWL91oqWMyYMSPuuvgrWvi3m/9xf4WBeNUOYtm8ebMt31bc96UqRWxyU7LVIR555JGEVRRifxeSpGKfJwGPsmLuPtsj9vvHJm2m8puXtjpEpr4XxCY9xv7Gifjss88KmjZtWuxvSaUKoEzhEUcckfB1rJ+/VJk/Mc6/romS0EjOi/e+jv79+xe7nv7vmyjJLtG+ToJivPf0r7cQQUCRYCECDpFCokhEjpjGJLKFF++AAw6wES86ouEjjAd+UZJrSHhzrYGJ8hCxjTdlSfST5w488EA7Td+oUSPb+YnpXQfTs0Qb8ViSEMW68N4uenrvvfcWShorDUz303qWpB/8tkTDiYgxPUukjMhW7FQtXdWIdJO0g/WChdu0UvavPxCx4zuwbVlvpnhr1aplI3VEwrLVLQ747Oeff94mSjGlzrrxGEludNoiUk00OLZbHNuGbnIXXnih3Ub8XkQ5+e38UX5mD4hwn3HGGfa93RT7/vvvX6TLHpF3ovtEqNl+/J50/yIKHTQy+b1IVOvWrZs93lKZ2uez2TfpOMc+RESU35L9mcj04MGDzdFHH21fyz49ffp020qZ44j14jfEW8yMAL8tyXKlhUgu0W72jXh2HdaFCDD7EduF1/A/3nuSHJlRKS1YjmhQw/o7u48QQWQblHCuV0IIkT5IsGratKm9zYWYi31x8HyrVq3sbYQCdgghhBAi6igSLIQQQggh8g6JYCGEEEIIkXfIDiGEEEIIIfIORYKFEEIIIUTeIREshBBCCCHyDolgIYQQQgiRd0gECyGEEEKIvKN0De/zFNpO0sKVgvTqjS6EEEIIETxogbFhwwbToEED2y49ERLBKYAAbty4cTp+HyGEEEIIkUF+/PFH2/k0ERLBKUAE2G3UqlWrmnTx119/md9++822SY3X3jJf0XbRttE+o+NJ5xmdg4OGrk3B3zbr16+3QUun2xIhxZUCzgKBAE63CCZcn+udJmhou2jbaJ/R8aTzjM7BQUPXpvBsm5Ksq0qME0IIIYQQeYdEsBBCCCGEyDskgoUQQgghRN4hESyEEEIIIfIOiWAhhBBCCJF3SAQLIYQQQoi8QyJYCCGEEELkHaUq4vbnn3+ajz/+2Cxfvtzer1evnjn44IPNdtttl+71E0IIIYQQIrci+JVXXjEPP/ywmT59utm8eXOh5ypVqmRatWplLr/8cnPSSSelez2FEEIIIYTIrgieNWuW6dOnj5k3b569X1BQEDc6/Prrr9vlwAMPNEOHDjVHH310+tZUCCGEEEIU4Y8//jCrV682v/76q/2fpXLlyuaEE04IROe2oJLUljn22GNt6znEb82aNc1RRx1l7Q+1a9e2z69atcraI95991274T/99FNz/PHH2/Z5QgghhBCibPz222/ml19+8USuX/DyXDzKly9v2rZtq02fgKSHBx06dDCXXXaZad++vd2o8di6dauZMmWKtUy89tpryb61EEIIIYSIAwHIUaNGmfnz56e8fT788EPTpk0bU66c6iCUWgS///775tBDDy3xdYjjTp062eWjjz5K5q2FEEIIIUQCvvzyy2IF8E477WRn6WvUqGH/Z3nvvffM4sWLzZo1a+z/u+++u7ZvaUVwMgI4lkMOOSTlvxFCCCGEEP8PxQgczZs3Nw0bNrRCt1atWlb4VqxYscjmwgeM+IUPPvhAIjgBcksLIYQQQgSQb7/91nz33Xf2dp06dczZZ59tc7RKYq+99jI77rij9Qp//vnn5vfffzc77LBDFtY4XCRlEmndunVSC1mIZ511lvUEk6mYDI888og54IADTNWqVe1yxBFHFPITr1ixwlxwwQWmQYMG9gfEk/z1118XqUxx1VVX2VER2ZBdunQxS5cuLfQapgS6d+9upw1YuL127drktpIQQgghRA6jwJShTUYAu0gwUWOgSMEnn3ySsXWMfCR4xowZSW94ePHFF81TTz1l3nnnnRJHHo0aNTJ33nmnadasmb0/YsQIc/LJJ9sfbJ999jGnnHKKDfVPnDjRiuT77rvPiu0vvvjCCl645pprzMsvv2zGjBljpwj+8Y9/2FrF+JJdEt+5555rhTGJe9CrVy8rhPk7IYQQQogg8dNPP1k/MFSvXt0cdNBBKf19ixYtrA5zlggqe4lS2iHi1QYuDsqk3X333WbAgAHFvq5z586F7g8aNMhGh+fOnWvFL/8Tyt93333t80SZmRJ47rnnzMUXX2zWrVtnBfczzzxjxTGQRdm4cWPz5ptvmhNPPNEsXLjQil/e67DDDrOveeKJJ2zUmR1szz33TOm7CSGEEEJkkrfeesu7TdnZRJW5EkE3X7TQjz/+aJYtW2ZFNX5ikaIIHjZsWDIvM3///beN0Do7xPjx40sUwbEl1l544QWzceNGK1CxOYC/HTM7Ad3paOCBCCbau2XLFtOuXTvvNVgn9ttvPzN79mwrgufMmWMtEE4Aw+GHH24f4zUSwUIIIYQICitXrrQBQKhSpUqpChRAy5YtrQh25dIkgkshgnv06GFSAW9u//79vczEkqD0B6IX4YyRG/GMFQJxu8suu9j3euyxx6z9ATvEzz//bJYvX27/ltuIYqYK/NStW9c+515D9DgWHnOviQci3AlxWL9+veevSWcjEMS/W4S2i/YZHUuZQOcZbRftM+E5logCE1gEbAxYUkujO5hFnzBhgv1bgoYEDONVk4jaeSbZbZWUCP7hhx+sAMVvWxwjR460/9MumWS3ZIszE4nFPkGi2rhx46zonjlzphXC3L/oootsGRCiwFgeaNyRjH3D72OO52mOfU0sQ4YMMQMHDizyONmW6Sw8zc5C5iakOt0RZbRdtG20z+h40nlG5+B8uzahhejPgAhmJhwhm6gjXDLsscceZt68eTawiBBmpjzq1+3fktxeSYngJk2a2Gjr/fffb3r37p3wdVRxQByiwBG1ycJ7u8Q4Qv4YuB944AEb/aXeMO+F93fz5s22VTO2Bjc1gOeFx6n+4I8GM5Vw5JFHeq+hykQstHsmYpwIItD9+vUrFAnGX0O0miVduBET7ykRrO2ifUbHUibQeUbbRftMOI4losC8Lwte4JICkCVBJBmrKvA/dtCon2f+/m8UPW2JcQjNK664wiaXPfroo4V8umVJoEv0Hn4bAuDfBcqj4Wu544477H1EMqH9qVOnmjPPPNM+hlUCLw2JeYDVAhHNyAp/DNBNhcecUI7Htttua5d4pUdY0onb4dP9vmFH20XbRvuMjiedZ3QOzpdr04YNG2y0lllq9Mexxx5b5s+gWxxC+tdff7U2VQJ6zK5H+bqd7GcnPafPD4I4pQoDotIVby4rN910ky3h8f3331tv8M0332xLsnXr1s0+T6Ic9ykYTZm0tm3b2rJpLhEOcYxdgrJo06ZNs6XVzjvvPLP//vt71SL23ntvW1/4kksusSKehduUUVNSnBBCCCGCAEn/2BaAoF06Glyg3yiXBug4RLb4P1IythL9ZQPiLSEC++qrr5qygk2Ber2I0TZt2tgILeXMELsuqsvzdD/p06ePvU15ND///ve/rTAmEkzYn52G+r/+UPzo0aOtMEY8s+BZRtALIYQQQYYp5rJ4QkU42LRpk61m5SKZRIHTBZrN5UAxm56OWfsosE1BElsCny8bjyjraaedZqOy7nF8s86a4F6X66zATMEUApFnbBQ07kgXeKg5weGhkR1C20X7jI6lTKDzTDi3C9PjTz75pA0IderUyRx33HFZ++ygb5tckantwmz266+/bm/j20VvpRN6KrjmGzQMa/bfXKwo7jPJ6rWUIsFETwmj05bYGY8HDx5sa/HiNRFCCCFEeiDLnsZOriTo5MmTvQQnES3Ig3Ld3QgmkhCXbpwlAihAIFIUwYCypuYc5cOc3YDObAcffLC2pxBCCJGmqXEEsL+WPRO3Y8aMsdWPRLQgcd+VFmvevHlGEtcoO+s8xuRgbdq0yeQ7pS52e8MNN9iKDDSc4MB0HUmEEEIIUbao4NNPP23b3ALTueTFAE2lRowYIQETIbAQvP322979TESBAXuCC1imWso2qpSp4wM/1Mcff2ybY8hkLYQQQpQNKgMMHz7cLFmyxN6nURXVjKiYVL9+fa/G/fPPP6/rbkRAR+FdddFaehtkCn/75Q9kiSibCAYOyunTp5u+ffuW9a2EEEKIvIXoHFWLqOUK22+/vbn44ottUydqxtJNlccAbzCzsSLckFuFhnK0bt06o5/XoEED06hRI3t76dKlZtmyZSafSUoEUxPYVYSIB97ge++911aPoNOJEEIIIVITQ88++6xZtGiRvY/opQZ+w4YNvdfgE6UOvit1RT4OjaFEePnss8/M6tWr7W2qNey8884Z/0x/NPjDDz80+UxSIniXXXaxS0kceOCBWS3fIoQQQkRBAGNvcIKWLqgXXnhhXEFE96+OHTt69/k7f/KcCA/YSLMZBXaQeOfKl3388cd2BiJfSUoE04wCM34qpn7+RgghhBDFC6GXXnrJzqQC4uSCCy4wu+66a8K/oYnCQQcd5F1vR44cqUS5ELJw4UKv/B2Bxt122y0rn4ulhuZhQEUK1iNfSUoE05KYA5IWxwsWLEj4Op7jNU2bNjXXXnttOtdTCCGEiJwAnjRpki2P5RpOYXcg2lsc2CHOOOMM6++EX375xVopiCiLcEaBW7Vq5dlcsoG/ZvD7/93/8pGk2nngQ2K65a677rJLzZo1beOM2rVre5mqfl8LPy6vEUIIIURRuE6+9tpr5t1337X3EUDnnHOOrQ6QDJUqVTLnn3++efDBB83GjRttJ7A33njDtG/fXps7BJBn5SqAUGBg7733zurnE3WuXr26WbNmjfnqq6/M2rVrTbVq1Uy+kVQk+OuvvzZXXHGFPeg4cBl1MoIZO3asXbjNYzyHl4nXslGFEEIIEb9F7owZM7z7Xbt2tXk1qUCAitJpRJCBxPR58+Zpc4cAfxGBbEeBgc879L8Jcmg3ugHnI0mJYEYLjDZ/+OEH2yb5qKOOsmKXDceCh+nII480gwYNsq/htfyNEEIIIQoze/ZsG7V1nHbaaYUy9lOBigInnXSSd/+FF17I+7JXQYcmKAQX3UCGmfVccMghh3jim5rB+djvIaU6wdgfbrzxRtvfmnZ7tG5kIWlu1qxZpn///raDnBBCCCGKQiLblClTvPsI2MMPP7xMm4rAFIIGNm/ebBPlXAteETz83eFIcnSR/GxTo0YNO4iCX3/9tdhSuFGl1Fue0UOtWrXsku0wvhBCCBFGyJ9x1ZZoYYsIKitcg4kmuyYICJrRo0crUS6A8Ns4ywrdAEs7A5CJBLkP87BmcG6GH0IIIUQeMmfOnEIR3HSBRZFEuR133NHeZ7o9X32eQYZZc1fFAxspuVa5ZN999/W6EDJAY5Y/n5AIFkIIIbIAbWpZgE5wLnKbLsjup8KEA1EjggMCE++tG7QggnNNxYoVbfMM2LJli+0QnE9IBAshhBBZjgIfccQRGbES4vF0ienffPNN3kX2gv774wl3NgTsEEGgadOm3m3XvCNfkAgWQgghMgxi1HlBt9tuu5TLoSULwpopbti6datZtGhRRj5HpAZRVn9N6HR4wdNFvXr1vNvLli0z+YREsBBCCJFhPv74Y1u5AajksO2222bss5wIhs8//zxjnyNS+/03bNhgb9OymMoMQaF27drWFgGKBJcApVdY4kGNYBYhhBBC/B/UX507d663OQ477LCMT2+7BDk6yRGFFLn9/f1l0Y477rhA/RzlypUzdevWtbfp/OssG/lAypHgCy64wPTs2TPuc02aNDG77rprOtZLCCGEiAQkG61YscITqP7p50yJGtd+meizOrjmli+++MKsWrXK3kYjNW7c2ASN+vXre4L9559/NvlCqewQ8bqK4D1K9JwQQgiRr/ijwGVtjJEs++23n3dblojcMnPmTO/28ccfb4JIgwYNvNv5ZImokMyLKLPy6aefFnos1hLhDrJc17wTQgghgsJvv/1m5s+fb29TDQA/aDagSgQJeDTmIBJJoKp8+fJZ+Wzx/3z//fd2AWYA9txzz0Bunvr/jQSDRHAM48ePN7fffrt3n2jvhRdeWGQjkvG4yy67ZPJ3EkIIIUIDXbjcTCndwSpUSCr2VGb4nL322ssGsKhMQUvc3XffPSufLeJHgY855pjAdtitn6ciOGk7BMLXb3Vw92OXvn37ZmpdhRBCiNAmxGXLChHPErFgwYKsfrYw1gdMFB522mknrylFENl+++29+tKI4HyxtiY1JPV7WAYOHGhHMrfeeqv3GPdr1qxpi3/TC10IIYTId0hI+/XXX+3tPfbYw14nswlT70SE//rrL2tZPPnkkwMbiYwiVIRwYvLoo4/O2ixAWaLBa9assdUh2G+zvb/mgqR+Ecp5uJIeiGC47bbbMrtmQgghRIjJZRQYqEWM+CYauX79elvCVJbF7EBN4I8++sjexpvdsmVLE3QaNGjgRa6JBueDCE65OsTff//t+ZuEEEIIUZS1a9cWmgp3JcuyjRpn5IbZs2fbCLyrC43dIOjU9/mC86VzXIXSZrtOnjzZZjySeRqL3yohhBCiMEyRalo62rz//vveVHiLFi1s7d5cgPjmswlg4Qvu2LGj9r0Mg50AEQxU5MAKEQbqSwQn1/qvQ4cO5pdffkn4GolgIYQoCpGhCRMm2GnSVq1amXbt2mkzRRAE5wcffGBvI0BzORVOWTYaNHzzzTf2uk0jBL/YEZkZAFGRAw466CA7ExAGatasaS00iPh8qRCR8tD0mmuusRmPiapDCCGEKAoXxaefftpeILGUvfnmm2bhwoXaVBEEG8S6devs7b333ttUq1Ytp+sjS0T24NieNWtWYFskF8c222zjdTMkQc4J+SiTsh2CmoNsKH7Y008/3Y4yNa0nhBDF+0OHDRtWJLrywgsv2LKSVapU0eaLELlOiIsngidOnGhvUyWibdu2uV6lyEJzMQQkUKc50y2y0039+vXNkiVL7G3OV8wiRJmURfCOO+5oNm7caMaNG+fVlBNCCBEfpp+JACOEYYcddjB169Y13333nc2vQAjTfEjBhGiwevVqWxoNatSoYasz5Boi0Y0bNzY//vijFTaUv2LdRHphNtzfHCNMUeBE7ZOjLoJTtkOcccYZ9n9O4EIIIRKzePFi8/DDD3sCGOFxxRVXmO7du3vR30WLFpk5c+ZoM0Y0ChyUwY2/XTPRYJF+8F27qgqNGjUKpYCsn2ed41IWwUyjMKo85ZRTzEMPPWSmTZtmC0L7FyGEyHewjj311FNeBR0uigjg2rVr2xm1M88803vtq6++aiPGIvyJjy4hjqoAtEkOCvIFZ54ZM2YUigIHZQCUCvXq1fPWWyI4DnScIarx008/mauvvtpmN5Pl7JbWrVuntMEfeeQRc8ABB5iqVavaha5zr732mvc804VXXnmlvYBQZ48kA/4mNtpy6qmn2osL78HFZcWKFYVeg0eH6AtZmizcdtEZIYRI95Tos88+69UJ5bx16aWXFvL+0s3rqKOOsre3bNlinnvuOe/1IpzMnz/f/P77717klcFOUOD66PypeD5pniHSB5ro66+/9mZ8/JH3MLHtttt6TTIYmFPpJMqUqnBhosoQpakQgbi98847zYcffmgXRDRC2/U5J2lkypQpZtSoUTaTmvtXXXWVZ/LHn4wQZ+Ty1ltvmXfffdds3rzZdO7cudCPd+6559rIDO/Fwm2EsBBCpAvOOZMmTbKRXQeF8nv06GEvLrFQs9UJE6Iu/gCACB9BS4hLFA3mOu0aeYj04J8FP/bYY3NWFzqdlogtW7YUWw43LxPjOJmnE8Sqn0GDBtlILycTDli8cnzm8ccfb5/v1auXeeyxx6xgRiwjemna8cknn9goMJCFzUgMUXzCCSdY8Yzw5T25IMETTzxho85ffvmljcgIIURZcNFcv9/yxBNPtAP7RNOiFStWNOecc4558MEHbRT4nXfesRnlu+++u36MkEHUzOXKkPjYtGlTEzT2228/a2F0UesgCvWwJkMSWHOJr0GywZRWBM+fP98bnNepU8dElZRFMAIzk/X1yJQmuotABTqtEFnp2bOnzVrEc0Pm7QMPPGCfp6gzFxh/lIU+3YzCqNWHCEZIY4FwAhg4+HmMri6JRDDvzeJw00dcrNI5bcn3dovQdtE+E75jiSnwkSNHeqWFOP9g0eJiWNK6ME3NbJaLHiOksZpRfjKd6DyT2e3CtcTfIS6I53PEDDk92ANJ4tqwYUOx7Xy1zyS3XRhYuJnnI4880h7/YbY21alTx9uXly5dWshPXhJB2WeS3f6lapvsoMwKI6CyRi0YcSB6SSDBQzV+/Hivz/rQoUPNJZdcYm0TFSpUsDvXk08+6bUhRMxysbjhhhvM4MGD7Q/HbXZIZ+pmhB5vJMNjxSWjDBkyxAwcOLDI4/iU0znVwc7ifGQkUwhtF+0z4TmWaIqAAOZcCJUqVTJdu3Y1zZo1s+eKZDjwwAPteZD8Bs6rY8aMsbkN6Uys0Xkmc9sFC957771nZwOI7lMWLdnfPtvstttuXjUSOsCy7yVC+0zJ24XfmVlmNAfBOLzAQf3tk6Vq1ap2XwZm2lP5PkHZZ5Jd51KJYA6cyy67zFoSOEmjuE877TSbaIZw9Edck4FILFMJ/D31h7E/kFiCEEYEs4MRDd5ll12s7+byyy+34XqivERRiB6zPrwWccr04sEHH1zoB4h3MUEwF3eR6d+/v+nXr1+hSDC1FhHq6Ux4cCMm3lMiWNtF+0y4jqUXX3zRnhsQP3wuNX/9tTaTpVu3bub++++3XZqI0lE6jYhiutB5JnPbhYoQiCD2AaL/tWrVMkHlkEMOsdduwL7hkjPjoX2m5O2ChYn9hgUvMJok7FSuXNkKYc5FzBqkoneCss8km9CXsgjm5Iw/F8uCPwkOcUiyGoI0VRFM5ISoCXAC4YSC3YELwk033WQjw506dbLPU0kCwXzPPfdYEQxMJRJBwcBNtJjpHpJNnCeL27HVIoD2z3i3EsGoLl4yC5/Bkk7cQZTu9w072i7aNkHeZ8gIx57FYBp7FYPx0jYh4O/OOussM2LECHsfewTnxXReVHU8pX+7cB0kUOMCKojKIJ/HqV1LlRIiZVQzQCxwDU6E9pnE2wWRiF7ht2cARFm0IP/2qdCgQQPz7bff2gE+Mx14nZMlCPtMsp+d8pz+HXfcYQ8efnA/RF9ju6WUFt4HLy7heJZY6wEbN57KZ/SNACYhbuXKlaZLly72cawWTFm+//773muZuuIx/DtCCFHWuqCUiCxrFy68dy6IwIVHZdOCD4EhZ71jthLrXpDheuo8nlxfXXc7kTok5jvbAMdtkEripbNpxrL/NgCJIimLYAzgjHpiS/m4mni0ZUwFIr1MJ+A7wRN388032wsLU4OE4xlZXXfddfYxpm6GDx9u/XcknfiT9RiJEw2mlBp+PEqpuYQ3anS2b9/eeot5HQu3TzrpJFWGEEKUCmaePvvsM3ubi1+6rAucl1z0l6SUN998U79QgOH65WA6PAxQJcKh7nGlgyiw81Y7K0SUaBDTPjmqpCyCibCCS0xzuKkg/COpgE2Ber0I1jZt2tgILeXM6EwHJIhwcUEU4xGmpjBl1Hr37u29B2XO6GCH2L399tutkMYu4Wf06NFWqGOdYMFW8cwzz6T69YUQwsLA3FnCOB/Gzo6VFixYzKy5GbDp06fbAb4IHly/8G5D9erVU8qiz3VyHFWUgHrBYa5kkCuwQbjqUfismYWOEvXzpH1yyoYNorMI3diqCpyo3YkgFWgrWhz4eUsqy4YwZikOpimJEgshRFnBSvXRRx/Z24gJV9IxXTClzuzV5MmTrdB+/vnn7exWceWsRPahDKcDL3BYGiTglyRoRH19qjLh/aSihUgOxK+LAhMAdH0MokTdunXtd+P8IzuED0Y84I/E3n333eb888+3Gyyd2cxCCBFEqFLjsqDJK8iEOGV6lSQmoHKOCzSIYEByOJWS3ECoZcuWJkzIElG2KDB2CDjooIMCXQ2ktFSsWNErLYsDINd1fzNFysNWypMxMsCy4CwQlBJzNgieF0KIqEINTJdky4Ui1hqWLogqUi3CnWexfYngQG6JS4oi+OPsBWGByK+z8CxYsCDpklL5juvs6E+IjSr1/2uJ4DtTTSuKpCyCaVV8/fXXWyHsX5wY7tChQybWUwghAjMF7ryARP8ymRGOvcxdiLCgueiTyC2IAjrEAYOU4mrtBhW8567RFZ3jfvjhh1yvUiigxrLrHkueEpbNqFI/DypElKqIG/5bKjBMmDDBJgbgHSExzVklhBAiiuCfpCySi9Qec8wxGf9MLBFcgAg2UCHHddMUuYNa9QhHZysoa2m8XEGyOIlxrkpEkyZNcr1KgYZoeWxZxChTPw+S4yqUxgvn/GoSvUKIfAIbhIvGNm/ePCviBxHsErAkgnMPg5EwlkWLB8lxDOYQd4hgmlKls1V31CCRkLbmrsJG0GtCl5X6eSCCU7ZDkAXZunXr+G9WrlxkuqUIIUTsFLgLAmQzI9wfnUMEi9xCuTonCHbeeWfbICOs0AUMMQeIu3idVcX/D378UeBszALlmqpVq9oWyiAR7MPfLjn2sXjPCSFE2KEkmvMCUg+2uJbr6QTPsfssmmc4P7LIDf4ocBSEkGsqBUuWLMnpugQZIuVukMDAJ8yDn2TZZpttvGgw9h9nAYoSSYVtMczT0S32ROAXvK5zkiLBQoiokWsvYNOmTe0FmPXgfOwSmkR2IUN+4cKFXtKi65QaZho3buzdZt9ybbvF/4PWoVuu//jPF9tI/fr1bWtwFw2uUqWKyTsRTLMKOrH5d4h4U4HsFFH3yAgh8g8G+atXr7a3EaB+4ZAtEUxJLmeJkAjOfRQ4TM0xiqNhw4aeL/jHH3/M9eoEEsoTuuoIbC/Ky1EnOh9oENM+OWpNVZI+gv2l0Pz3Y5cePXpkal2FECLrcF7zN6rIRUY4IthBdy+Rm/rQrksg5cWi0hiqUqVKXpkvZhuogCIK4z/+yYnKlyhwPiTHJRUJpiOKE7cjRoywOwAd4hzcr1mzpu2cdOqpp2ZubYUQIsssWrSoUCKUSyTKJtWqVbOVKEheYsqaJD1Zz3LbHCNKLazZr10ZPnznzZo1y/UqBQYGnS4hFW8+JfGi2j0tHnSNK1++vP3OUawVXCHZBhksTgQ7i4QQQuRbFDhXUSCiwYhgBDDT1v7osMhuc4xMdQnMFdh7nN2GQVbYRfBvv/1mB65uQby52zznfsdkFlqWIwIZhHbu3DmvosDAYBshzLajfXLUBuApfxOV6BFC5Auc71xSMFPGuWxUQb1gNx3PekkEZ4958+YVqgwS1uYYifBXOghD5ziEGKXqaPfMseAXuSyZqmLw5ptv2uOQfQBvLL0SDjzwQLPTTjuZqFsili9fbn3jWGbwReetCOZgITpC0XguDvH8Q36rhBBChBV/FJhk4FxGgfyiV8GI7BGl5hiJqF27ttluu+3s9ZxZBr5zUCKe69atsxU5ELx0t+P/r776KielAtkuiG+WWF2ETYJqIYhiGulEqYpC/Zj2yXktgvHHdOnSxSsTE0usX1gIIcLITz/9ZLPCXTksciNyCXkXFK8nIkk9V6IyUahOEHS45jkvJLaBKNaH5brNd/v6669tFBULAPt8tmE7E3X3C148ykGH45Hl1Vdf9bbnXnvtZQ4++GAbLWYhlyAoA4uyVoiIEimL4CuvvNLrNS6EEPkSBc614OQCSjQYkUDEDsGgkpTZb44RViGTTHIcIhiIBmdaBBPJnT9/vrX4fPzxx/b/dAssbCtEMf0LSaaxDb7iLawPApxBJ0lhDIqxYSQDf0+gkGX06NH2MT7XL4qJFtMIJwzUj3CFiJRF8Jw5c7xRTseOHW1LvaieFIQQ+dsUgQs0MK156KGHmiDgRLCLUEoEZ685BiImCs0xEuGvfU1U84ADDkjr+zNoQ1i6hQ5smzdvLnOEcu+997br7kQuj/E/lRxKW8EDoUeDCLy+CNUbbrjBClsqxXBeoG44xyGDhmS/A9H1t956yy7AoJpufQhixHHXrl1zPtBOBNuA8yCzBGybINllsi6CyZJ0BnH/6EAIIaIC3eFcpIjoX8WKFU0QICnH7wuOoj81SLz77rvefkBzDHf9i2ok2FHWphmUksPK8OGHH3pLWSKIHH8kopGQRnKqWzIVrZ46dWqhWSDqQgOWKBYiwlSZoMYyg1EEPeIYYcxtV0qvOLAzuWgxeuqss84yQaZBgwbWHka9bHzaLqKedyK4U6dOZtSoUXnTLUUIkV8QsWF6FkgWOvzww01QcNGtTZs2WREcpYhM0OBij3gDRFDLli1NlCHa52pRM/WfSi3cNWvW2GPmgw8+sNvsk08+KXXTDaKvRNwRuU70UrINwZkN+O4IWSD6ecQRRyR8LetEQhzL2WefbR/j2OTv/VFvKioUB9HgoFO/fn0vR4IBTd6K4Msuu8y88sor5vTTTzeDBw+2tojYKIl/RCmEEGGCyjdOABD9QwgHBecLJi8DkcbF1XX7Eunlvffe86a6scNEqTlGIrh2I4KJZP7888920BULxwYDMMSdi/JSraG0+7MrNeZPIMulLcAfBaY7XKqzQOwnNFNxHQUZqCKsnSBmsBAbLQ6LCHYggrGh5KUIpiucKyBNlYhYeC5Z87gQQgQJLlhEsdy5LEhRYL8lwiUnI0YkgtNP1JtjJAJv7aeffmpvUwaMaW8EG/YI7uOB5f/SliejuklschiPBQW+pzu2iHSmI/rP/oN3n8U1HYuNFhcXbQ4K9SOaHFfqth/OJyWEEFGBi+Dq1avtbSJSQSyCH1svOAwX0LCBEEQAAtPxlKeLKkR2meYmmoud4e2337bf/bnnnrO+1bLAMUREFMFLNB1bQ1CTv+CNN94oUxS4tNHiMFCnTh3bKY4BYpTaJ6csgknEkAdNCBFFXBQYiFIFEQrV41ElGidfcPpBAGL5c5AYGSWweFDdALsHrZIRvunosIZtiEYRftEbps56NP/y1wUPSkWYoFCuXDlrj8Ha8csvv9j9KFs+7UCJYLKmhRAiahD1cuXHiHiQ7BLUixENG4jcIdjwcEY5UplNmOEcO3as9VsDCVphb0/N1DuDO0QvJU6Zfi9t0lpsZBChiOjlf46XMIsifxT4hBNOsOcAUdQSgQjmOCEfwV9WL6yU6VdmJMAJmNFemHd+IYTA70jZIyDpI8iJUAgzl4xEiSaJ4PQwa9Ysr2EEVhgSwMM280nCFWKX8m4IX6wdyZTsKq6lMslr1LTFzsBtlijtc/icqQsMfC98y6Lk9sl5K4KZguvTp48dOeEPYcR04oknmn//+9/WAySEEGHDJQQF2QqRyBccJm9hUKEawmuvvebdP/PMM80OO+xgwgAiF9H78ssvmylTptjE9dLYbBC6JKqR+MQggEQu6uRSPi2qkVGimv6KEESBo1wPuiw0iGD75JT3akLgVIhYuXKllxzHAUjPbEqlUP5D2cpCiDDBrJarDYq3ETEQZIjAuCQVIsGibHANe/75573KRviAd99991BUsED4Tp48OWXhS3nTww47zFZA4X9XDg2/5913321vuyTRKEME2B1DRL2DPgDOJQ0kgo2tDewKPzNNxKiR3trOIzJkyBDzwAMP5Pq3EkKIpKFrkyv7hA80KB3iEsH6UdOVizeWNARQVIrX5wJa2RIJ5prGlG/79u1NkIUviXsIX5pUJOsjZ792opeZg0RJa9gBiIDji/7hhx8iXQmK7+b3Ardt2zbQ1StyDRYxzjOcb6LSPjnlSDDTRXzp8847z4pdt0Guvvpq88wzz9gDUyJYCBEmwlAVIp4lwkWwyGynnasoXSQQDy0DC6LrdP4K0iCIwRneXmZbub4y6EkG9gfqGyN6SVyj+1kycH1npsG1yOX6nuzfhg2+45IlS+xtZrCpbiFKjgazT5BcySAsTBVA0iKCXU/x+++/34s88D/3EcFLly5N/1oKIUSG4EK/aNEiexsfJM0owgDrOW3aNM8XLBFcut/+hRde8O4TAfYn/+QKrqNEp1nw+rpqFSXBAO6kk06yC80ZSosTwW5dopAAlUwUOOxRzWxQv359r6EIyXF5J4IZKeOfYwRALT2HG51G1TwvhIgm1Ex1TQGIBIVlOhQ7BOvKuiOCReoiaNy4cdbOB1Q+yFVNYK6pRHunT59uha+rUJEMDH6c8E2XWGXfclASK4og5FzQDmEX1JKIQe8ct1/It1vKipWEEaYOyZwdMGCArVfJdMLtt9/uPS+EEGEhTFUh/NAwg4x+Zufws27cuNFUrlw516sVGqiXO3/+fM/reMYZZ2Q1Eoi4RPAifCnNxu+XLAcccIDp3LmzFb5+wZou/O8Zxdnd2Chwu3btFAXO0/bJKYtg/FJUgEAIuz7YDk4gZ511VjrXTwghMgYzWM5XS3a8P/s5LJYIZ1EjGhz2qEw2f/eJEyd69xGUmW6Rjbf3/ffft6KXplPObpAMXFuJ+GLXQPg2adIko+tKYlytWrXMqlWrrNBxpVCjApVgnIDDNkJrbJEcJE7SF4LZiyi0T055ryYBjpMHPqVYKJ12zTXXpGvdhBAi41YIByIjbJ5AkuNmzpxpb0sEJwf2kTFjxnjVQGjxmykRxCwpopeFaybd25IFuyE1elu1amX/z3ZzCqLBiOCtW7dawRj2znn+399fF1hR4NTAgkUSIZVDGEyyTwe5sVDaRTBZs2+++aYZOnSomTRpki2LRgSlS5cu5qqrrgpUVq0QQkStKoQfhAnCneld+YKTAwsC1TSApB4iq64+cFlBEFBpgkgvn5Pqb8JADNHLwu1cNm3AX4xlBJhtiIoIZuCLfQiwc8rCmTrYsBDBQDQ4zE3SypXWi3bdddeZd955x7bu5H/uU2Q+VR555BHrb6LeMMsRRxxRqGsPbUyvvPJKO2XBaIN2pvyNH3bo7t2729EJnjhaHr744ouFXkMiH69hyouF26XpqiOEiAZEt9yUKBfDMGY5c050zYnwmFK2SCQGMUcQBxg8nHPOOaW6bsW+59NPP226detm7ShcW5566qmkBDCVlbAVUlZ03rx5tgzatddea6PTue5a5vcFO8tN2FEUOD008lUeCXviZFKRYELeGPc5aXTq1KlI9jTTJdQvJBpBXcJULiZszDvvvNNm5sKIESPsSYEIzb777mv69u1rp5JGjRplfVCY2S+//HLr3XOeZE4669ats5FpfEzPPvus9SbTwc5Fd84991xr8KelJPTq1cv+Hd12hBD5HQUOc3kxInSucD3T74psxQf7w3PPPedVAmnTpo0d/KQaBeZ6x77DlDqC2pXXS8Xbi72BhetTrsVucQlQrBvd9KIigvndsHg4P73THSL1SHBUEieTigQ/8cQT5tRTTzUjR46MWz6IAwXxymuefPLJlFaAhISOHTuaPfbYwy6DBg2yfcrnzp1rn2dqqUePHvaEgQhGvFLGCIHr4DVYMVq2bGl37H/+8592hE0Cn+sGhfhl3Yg0s/Cd6LqTSnKCELFs2LDBdnByXRRFOEAwEnkDzmlhLpLvn6aWJSIxnO9pCeyinIjgZKFyAzOU/fr1s8KVAMxDDz2UlADGy3v66afb1zMVz3oQ7aWBRVAFMJAI58QO7ZOTrVUcVBi8uFkAUF3g0oMF1iVK5kUk2EVLL7vssoSv6d27t3nppZfsa6+//vpS76QULueEg1AFIstEeHv27Gmjv3itsGD4u9LxGvq+E6VG/I4dO9aO+hHOTiRjgaBlpIMuOjyGgEkUOeE9XPIEuHqSRA7S5SFz39stIjzbZcGCBbbOKD5AvPB9+vSxMxHZIOjbJlcku10Qi662OX42psTTeUxn27vpWtvS/SzR98jnfQYPsAuskNlOOTS2Gdsq0XbB64hoYuEaQjZ8MiBsseQdd9xx9hrEjKY/eBSm/QwRzD7lPOdhnmXAo+0GQRzzpZkFcOTzseQXwkSBV65caW2rzlYUlG2T7G+blAh20QUirYlwAtOVG0oFajUievGzEQUeP368l61LAt4ll1xibROMPDiZENFF+DoQwNgfGHHzGsq78B7OrI1nuE6dOkU+l8ecQT4eQ4YMMQMHDizyOD94Ogvqs7O4UXaQIwPZJqjbhYvh66+/7iWNAFOGTLVecMEFWakwENRtk2uS3S40JuA3Ay7sHNNhhXMRbW1duTfyH+IlKOfzPkNlBvd7Uw2AC7b7zd12QejRpMI1rGAGMZVKDjTaQPRyLSO/xRHmCCrXVCcmCD75p8HDBNFKZoOxwnB+5rcqyzGfz8eSf99w2pCBkivbF5Rtk+zvm5QIdqMnRtDFJcu5aZNU4SJEwXoS1YisYX+g7A9CGBHMCJ5oMCO3t99+23qC8SudcMIJ9u+xP3DiZ8ROJG7ChAmma9euNmFv//33t6+JJ0w46RUnWPr372+nv/yRYKIuCHWWdOFGTLxnvh5QYdkujHwZdHFMOKHBurGuRI4Y0FEqMB+3TRBIZrtwUeeCzu/HQhKSO3+FFc6hziLGuTBe6+d83Wf4vblI81tzDTvqqKO8axnCmOsLOS1cc1KZ2t19993tNQhbBZHfKG5TLIoElth2nPPSed3LFszmMkPN78PC71XWiHa+Hkt+CDI6SxnnHFejPCjbxnn/0yKCsQ0gbjnJJrrAuxOwfwScLJyQnEEdn9QHH3xg7Q7333+/uemmm2xUF6sDUEkCwXzPPffYE9DixYut14ri10w7Af4+BPB//vMf8+ijj9rs6XieTQzyhPQTwYUx3sWRk0K6C4e7AzRKBcmjtF04oLhIEgF20QT2W7yBRIEef/xx+zqSZTgZZKPSQFC2TdAoabsggJl14jfkt4pClzXOn25mguQ4xEs88nGfIe8DIcTvTVCE4AdRQY5lIr7JVglim2Gjw0vKQlAm6hBU4vhgsMAAgX0nbLW0CaAxS8J6E8Tit0uHOMvHY8kP+7/bF0jM9W+HIGybZD87qVdxocCLi5mf6gyxo0E8vJRIc68tK5ykOGlx4LHEWg/YuE7lu7B7ca9heorqEXTrcZYOpkN5LBtROxFuGOUS/fVbfbDnUF6pdu3a9j4XRyJK7Ld40i+99NLQXSzyhahUhfCj5Ljif2+ELt5F7A4EVvy5HsVBjknr1q2tcMLqUJogT5jhHMa5jmlvrrUEw7KV95AOvvjii0JecM7Z+Rq5zVRy3F9//RXqChFJiWCiXYhghCMRhiuuuMJGXTlASA4iEou3lvtUiEgFTkgdOnSwIzQy7enkw2cxUueEQ3IBApt6mIw8iMZRpeK+++6zf7/XXnvZKAiig+gwPhXsEETkyMIFagvTbhJv8WOPPWYfo8oERdLDbPQXmYfpHiw6rv4q+ziF7GOjCVQ4IeKEYEYsc+J1yZ0iOPA7cmEEcgeicvwzG4FgQ+wRCY5am9tUvYAIX2YnmVVMJanNVY7AN3ziiSfaoEm+bkeHE8FAg4SwiGD0hL9fAJWowrLuYaB8+fLWlkr5PKwyYe0cl9TRjXjEnsDJFbF76623FnreZSYTjbjoootSWgFsCtTrJZyO7QK7AwIYkQGIYry5FCJnSgMhTBk1qlEAXiX8XDfeeKPdyTkBIoop2YYwcYwePdpm73NyAzrcId6FSCSWaA/uT35DZJx99tlx/ZYk2pBxTuk9YJ9EYIWxAUOUYdDuEqSwTUUlKsTgjP2SspBu6jofpuu59hCFQuwiellIaEvWD+ggqIPoJSBD0ESzOP+PPxkOwYP/OQz7BcELlxxFflFxif2idBC8dDWkw9o5LikRTMQEQcBJAhHsRK8fRgS8htemAp11igM/77Bhw0pMUGCHLw7ECA03hCgJBntUenAltJxgOu2004od6bIfUiWFGROmW9knL774Yl1QAwT5BGFuk1wcBCFcbXRKgkVRBCNqqLVLpJeF71uaGt1Ed7HCcU3D7sCsIza/fI/6xoPSpA7XKjfocA52Mz78rgQoNLDJ7ADpp59+iq4IBhIKSD679957rdHcTY8QfcAuQRUFpuSECDPs11hmXCSJCC/7N9GPZE6iJHBii2BaGv8hPnR/fWqR2+lRkuKAc1XURGKsLxgrWZjB0sGx5AQvAxh+v1SjvP4LNtFAEqqxNDHz6D4nzCXyMg0Df3IfmPIm2sdMQ7wSfEGBhHdnhQQqRYWxqkUYaBiBznEpDXuJpmJFYBEiilAj1F1kqXuI/SEVSwOime5Qbobj1VdftT56DRBzDxFEN4tFFDhqkSGEChd7BB2+dPbjdNYzz4atAaHLgujl98JnWBqI6DIowO6CDxRrHKU3RemnvRHBrhRkUAeQrB8WSucBJ2EZe4vIDHUjkBynuR8h/gtJba6NNqIV33lpRARe4BYtWlifIt5ibBF45aMmusJGFKtC+GH/QvhRq5r9Duuafyo7aIIXketfki1VFg+OV8prsnDskVvyzDPP2PrAgOVBlE0Eu+MHS0RQRfC0adM8jyqDQldaVWQvOa5igGcJ4iERLMR/IanGRQq5mJYlikblEaZuKcPH/7w3F2eRGzhBOz8jJ21yDaKIE8HOEpFrERxP8LJ+DDhLC5EnontE8xnMcKxiy/MPMrG+UEMeqBhEhQNRNhHscCIzaOCDRwQD527KoYW9CU4Yk+N2CegAKRESwUL8txkGkVvgYlpWwYqPDlvE008/be/TsYjEOSpMiOwT5YQ4P/4LULbFCnVkmUmhOsOiRYvs/yQnlSXC60qW8Zs50Ust+pJKMVHa0A1o+RvNwpQNBo1E+PADBzE5jpkPbBDuN6e6lAY+ufEF7yIRLET4IFrrLtbUnk6HWOV9iFIRBeYk/dJLL5kLL7xQF+Qsw4Ux6lYIB1Fu59HLlFjBd0nUzQld9z9VVeJVDkoFjjt+HxYnfInkpopr5+oqu4iyT3sjdvjdqZqD7zxIyWauK5zL5SDxUeSmQkTYUCRYCGNsFQdHOutJUruaKhHYIhAL1B1GGIvswYmZjHFg2jzK0XgEMBEwxAoWELp5lrUtNJ3W2G9dHV6qBLnmMWWB3wHvLgvVh/ifqdWyRm0RQwjyqFtfsg0RefYrYIBF7d0gwP7IfgnYH0hmDktCaBSoG/LkOIlgkffgH3Q1JakXms5sYqZt6aI4fPhwL2KBLcKVZxLZHeBE2QrhQEg6sYIlghmJVKK8WBpc8wnErxOU6RS8RGcR65mwKSgKnDkR7GC/CoIIJrjg7wpHEyw1KMr+LEGDBg3swIiBdzoGyIEWwYywWFD9sZCBy0nNmdOFCANc7F1ZtLImxMWDi8UhhxxiBQUniPHjx9tyTfIpZh62t7NCECXKh6lxv1jhwlScCKapC40FnOil+URZauZy7JCcx2cymGShG1umBG9J/u8oW19ymRwXBF+w6wqHFx0YXGmWLTc0atTI2ydIjqtTp46JdCQ4ke9rxowZurCLUMG+7I8UZqqCA7YIfMcu6owwC0P70bDDdkboOUFEHed8FytYQ95++20zdepUM3PmTE9EpAr1dxG5fsHLLEdJSWuZhO5xy5cv9wYDigqmD6L5VapUsecwIsGcO3M5kEdsYTFzM3gkIiuwEIzkuDpRF8HxoPQNaCcUYYIySqtXr7a3mzVrVqoknGSgnThtl0eMGOE15ZAIzixcpOfOnevdp3B+PkDNXCdWEMHMcpC49sYbb9iFEmWpJrDhK2Z/JdLGrAYVGqjDGjT8UeB8iPpnE67tDLAYxDPDwmAql2LH+YCBLoCcY0UwkuMODlGAJykRPHDgQHP77bd79zmB4gOJR5hGAEL4o8CZbm/MtDDlY/BYErHCzyZvcOZgO7uoINvdf6KOulghIYwoGRckZjdonJEKbC/ErmtAQbQ30Tk/KHBdciKYbSARnH6Irrv8CQZYubreY8d0NidKt+m3zn1yXMX/ltALW4WIpCPBsZGDRJEERmRChAEy58ksBqIIiNRMQ7TZJRrRzQqhITJDvkWBmdF48803bbSXmQZnAykJBCNWEQaBLtIbxmAGF183q7PbbrvZKXKRWb95rjy4zGw4Gw/n7VxacITxOsexT3AMMlMQpBJ6ZRbBeIFcAWQu4Jw0/QcD95lGPvLII82AAQMyt7ZCpBES1VyCJxd+yrxkQwS7xFGJ4MxBcpezaDHAoSpBFGEfwtv7+uuvF+p4mIy94bjjjrNNBdq0aZMxG1A2yZda0LnEJTiyn+UyOY5zt0OdOIOXHLd8+XKbMxAGkrrqX3311XYBlzlPS04hwgoncdchLt21gYuDwaSbNkLA5Dq5JB8GOFwkw9bPvrgSZlRwQPQS8XWtgZO9SCF8O3ToYAMWUWopy3HkBj1EpfAsi/RDYilT39hrWDZv3mwqVaqU1U2N190lxBGgI+ovck9Dn92MpEWqdYSBlENfw4YNy8yaCJFFnC8XKOnEiT0bEG3m86gUgSeYuopBTDAKuyCi7Fe2vN6ZhmlfqjkgerE7uCn/ZCAag/gdMmSIrd6ABYhpymzMemSTb7/91h5PsOeeeypJKoMwC4wAJuGS81i2BxxE/F1JSxKw1BgjeCJ4+X9zMcJAymdC6ps6yA7dtGlTkdf4rRJC5FOHuGQtEVw8gG5yEsHphQj7mjVrbIR9jz32CM20nB/8vPh6J06caO0O8c6z8cD6cfzxx9vcDP7GJTFhf4jyjIMaZGQPInzu/MmMSzZFMANcf1UI5VQEh3r16tkZN2YHiASHhVKFA+666y5z7733xo1IcKKN10hDiKCAOHAXTab3sj1tgwj2CzampkX68NtcwpQQx3lz1qxZVvhOmTLFrF+/Pqm/YxYD0XviiSeao446yquFzHZwIpi6ri6vI2pgEXFWCC7CQehkFmWoBU3SIfsntoR0tOZOJfnRVTpp0qSJAggBoly5cjY5jllWWpeHJTmuVHaI/v37e+Z4IcIGZZTw5LrptGx72mgxSTYzYhxPJ1N7mtJLD2vXrrXRdab78QsGXRDx2xNVQ/i+8sor9uKRDFgbSGpD+JL0F2//ic3kjyr83q5SAL93lLzOQYR9jfbjNFphAML5lMFXNvBHgdUdLng0atTIq37EgAVrUuRE8NNPP23/5yLOiQcxTHF2Tt5cdFiECDJ+v2i2rRDuIkIyB+XZEMJMHXHyEGUHQekG5yTEBXFwwfoxE4HwnTRpUlI1fEn2wtuM6EX8JhPVpcwZUWEiMojgqAYt/FYIVYXIDtgQEMHOEpENERxbGzgsiVf53DRjzyiKYC7cCF8SNNw0Lsk9t912m3n88cfNa6+9lon1FCIt0NLR+ZXofkRUNhdgiXA1irFESASn5yLprBCI31wMcIpj5cqV5oUXXjBjxoyxiVwlwXkW4XvyySebTp06pVzGzHX4IlLK1DVL0BtepAozOu44QvCH4aIbFf8n507OpZxTSTLOdHIx1h7njUcAqzZw8EVwGChXmvqbwHSIS7RgSu+6666zB8I111yT/rUUIgIJcbG+OgciWJSdBQsWeOcnpsWD0I0PYU5VhwsvvNBO3w4ePLhEAcy5laAC++q4cePM+eefX+o6vlG3ROBJdU1BSNCKWtWLIONPSvPX7c0UskKEY3BU4b/HYFhEcMpnDHrSU4qGSAuZyFgiOFk7AzSJHUIEES6Wrq0qvsFcttqkYgEijWOJmtuIJV3Ay8acOXMCUxYNoUvEl8gvEeBkPL5dunSxCwk/6YJIsF8EU54vSvgbZDB4ENmD7f3qq6/aIBi/Q/v27TNmP2IW48svv7S3VRs4+Mlx3377rS2cQOQ+6BH7lEUwPjMu3HxBfI3z5883rVu39nb+bGWJCpEqZJDjjwSSiVwWfS5gFgVLBBEUpnRJJlDR99LDLJSLsBI1zcW25IRPchvi19+yORGIXawOLJmaxvdHgqkQESU4ll3TBIIwOn6yC9t8r732sjYFNAFJvv4ZrnSCyHaediLQUS73FwVLxLf/PRcTDfZXQ4qECMaLQ41TImqcvBEWRNhc56uOHTtmZk2FSKMVIteRQnAi2FkidBEvPX7Rie0gWxdJ16ns2WeftYludLMqDqIinTt3Nuecc45N3Mv0eiJUGBQQtOCCRDZ/lOwvrhwnszpBTIKMOghSV4aPc1kmRLBqA4fbF9wsaiJ44MCB5rzzzrORi2OPPda28GRKhJM5mcv33HNPZtZUiDJABr4r3cJ0jX+aOFfE1gvm+BGpwyDcDSawlGTD5kIzjvHjx5vnnnvOEwElTR0jfLE7YCnLJkSDEcEIRiLmQfBKpwNnbYJcWpvyGWw8rtwjs8KnnHJK2mfYXOId7LrrrqFsfpNPNAxZclzKIpiEE3/tzZdfftlOhXDxkRVChKUsWhCm0xAj2IvwjDJVzfRuLi0aUbC5IIbIVcgEeB/xHSN8J0+e7CVkJaJGjRrm9NNPt+I3l1ULEMHOO4ugoIte2CEBkqoXQInOqDYCCTpc9ylLx3GBrQshzAxHOlGHuHBRt25dL7+F803QqVCWCwL1GSmPRt1KIYIKJ2dmLFx9SRpkBAWiwYhgjid8VEFv7pCPCXHMIowdO9Z6fd1sQiIYXNG2+Oyzzzbt2rXLeiOWknzBYYjMJAORf44ZQIQFYVCbr3A+dccgv0s6RTDnbhfx51hSbeDgU65cOSuEua6hD4OeHFcqEUwUpFevXmb58uVem2RqBnOxoFYwLTyFCApEJ4JaXxIRPHv2bM8SIRGcGkTQXbSBaTjqLdPGNR0X37feestGfadNm+YJrkTwuQjfM888s9B0YBDA/kNkhu8UheQ4PKJ+f3+6I48i9UFW7dq1zapVq+xAnsZZzIJkojawZsrCQYMGDbyqOEH3BaecSUBE7dRTT7UCmJORy9ikQsT3339vfXJCBImg1AaOB8lwLopFdrUofRT4iCOOKHNEkIHIoEGD7H7Ss2dPM3Xq1IQCmMgUHl+EMuvRt2/fwAlgQAC79UKguBbDYYWSgggud/zII5pbOOb8s2tu1i0dyAoRThr4mlAFffYpZRFMsXciCrEnnq5du9r/33333fStnRBlxEUngGhF0OqkEpV2AoWBpWv2IEqGCJFrmUuEqLTJUdQgHTVqlK3acNxxx5mHH3642Nq++HtJEGbq95FHHrEJwkGvTOBPBA17NDhoVV7E/1ki3AAUEZyOFt3kGlGJyvm+VT0nPNSvX9+7HXRfcMpn7nfeecfu7LHtkV3iR9BVv8gvnNUgSAlxJVWJEMlHiRiQu1JNNEBJFqK7nMuuuuoqW7nhhhtuKDaCRdLvueeeaxOBsUdcfPHFaZvyzQb+xLEwi2BXhQBIgNx3331zvUoiRqTiAy3JO58Mqg0cXurUqROaznEpe4LXrl3rNRvws3nzZi+qIkRQLphuOo2EuKB6B6mtOWPGDE8Ek+gjiodIk7828OGHH57UJqNrGkludHJLJkKBuEb8EiUOc/Wb2M5xYYWBihv4EH3kuBbBgN/DDeL5ncrS+TC2NjC1v0V4KFeunLVEMOAOenJcyiKY6AdThW6K2UGheDflLEQQ4CTqylghZjJVOquscLFg1EyCqSLByYF/2vlCqR1KNnIi8MCSzEt1B7+HOBGcwyhtdtZZZ0WinJiL1NE4g/rGiH/X3CjMCXFB8/fnOySuTZgwwQbEsCkxcCztIIWBmrMkcXyHadZF/B9OBAc9OS5lO4SLuHTr1s177PLLL7fVIjipHnXUUSm9H546ospVq1a1C8ktfqsFHskrr7zSZl8zkqA4N3/jIBmPz423EO1xcPLv3r27rc3Kwm0X1RbRgwum3wpB9ZKgwoXCTVeTuMQiUkuIi2d3mDVrlrnmmmusV/jqq68uVgAzCOnQoYMZNmyY+eCDD8wtt9wSGQEMnA9dNJiojBtAhAnEO7554HipV69erldJ+MCO5EqYsY8tXLiw1NvHNb9xAQwRPhqGpGlGyiKYDGg33eEiCY899piNuHGfi04qIG7vvPNOG7VjocoE7Zhpiek+b8qUKTZxhYOK+/j4XOSZEzsnRv9C0gpTl1zUHExpUm+Q92LhNkJYRBP2FbpkObtB0C+Y/lGyawIg4kNrYnd+oPua3xdKJP3++++3g3EiuQyEi6uGwKB6wIAB9nz25JNP2tq+UZ1iD7slQlHg4OMXrH4hW9rawAjrWOulCJ8IXhrg5LiU7RBkQj/00EOmX79+hTomsbNy8YkXlSkOpkz8UJ6ISC9+Py5uRG969OhhC9ADEWdEN4IZsVy+fPkiAocybVwAmf5zggjhy3u6bOInnnjCruuXX36Z025OIjP4q5SkOjuRKxH8+uuve0JOWe/Fd4hzZcvwCpKHwKD4xRdfLNRKNxHVqlWzZR45R+y3336hswWko2kGIjhMPkuuNe63pRKIhFEwwbrA8cUsK9dWBqyptglngOs6QBJZTiXhVQSHunXr2oACg5ogR4JL1SzjsssuswIU2wI9vfmyRF39teFKw9atW23khmL3TkwfffTRZtKkSbZmJ+9PAhFlUx544IG478Hok5Plf/7zH+8xhDQWCL+wwNbBY0yZSwRHC/ZJF03FS7bXXnuZMETpuLhz8sfvGkbPZrbAb8i5YtmyZXaGiPORS5YqLlGDgTSlHIn25mPRfSIzbp8KWySY39wFXbC3SBgFE44zEuRoNMNAlQoPBM5SQbWBo5Uct2TJkkAnx5W6bTJf7qKLLkrLSlDyBtGLACB6SyTXdc4aOnSoueSSS6xtAt8eG5ZpS8RxPJ566ik7xen3gNLJjpIdsfAYzyWCk64/2u0qX5DAxJIuuKC7RZR9u7z99ttenUoGO5yMS+r4FQTwOS5atMhGT5g+8tdajCVf9xkGuUR8yQXwH5uJ4Dxy2mmn2UG7/xyQzuM3LBCVob47+RHYxgg2hEVMMovnjmmm3HX+zQ6lOc8QpaeMIOCvTyUfw9UG5rcmmZPZiyAeq/l6/k112zBLz7kaEMPZTI5Ldr+pkKyoSIVUR35EYoneMoUybtw4a3+YOXOmvYAhgjkBEg1GJLAuJOIhEGLbMzPSePbZZ21SSyzxomolRduGDBli/cWxkKyXzuL47CzOt4i9Q5R+u7AP4B0kMkhHL6LAYWlAQaTO1UD9/PPPi51GzKd9ht8UO1Oydgei/yeddJIVvv5ZgLDsB5mCfYbKF2Tdu0ok/vrBQZ7ZcdWIuKgyg5fO3zKfjqVsbBuq8DDgZAqc6gDsZ8nkZHA9JoLsyq1y/U9HC/RMoH0muW3DudjN0jHDmc3cnGTPEUmJYKYRk52a5XWpjtwQK26EgE+N0SN2BzzGN910k40Md+rUyRtlciG85557iohgLpJs/PPPP7/Q42x4TqSxkCFdXGml/v37W++zPxLMtDXRauc3TgduNMl76iRctu3iGh4Q9SIKHKaWqvjfXASFC0hx+1g+7DN4A2lJTNklouPFQUTzxBNPNKeccoqts4wvMarbpbSwz+DZJNLGeZrE0TA0m0AYuWRF/P2pekxLIh+OpWxvG869BK4Ab3BxEUBm6Rj8Y3VkZtb/W6fzOptOtM8kt2343d3vSdWjbP6eyc7+pmSHSEcrxGQ/h6lORhAssVFXDsZ4XxArRJcuXYrUKsZqwTQLEUJXW/K9996zjxU3VcOFNd50IbYM1w0lXfCdWNL9vmEnle3CPsHv6gZsxxxzTKi2JxYjygQi+Jg64nsUd+GJ4j7D6B3Ry4yOa4lcHPj8zzjjDBv5ZdsxAOc9orZd0gXTy66EJJaboG8jzv/sB6wvF1OCJJlY5ygeS7ncNlhWqM2NICKRlQT42Os4xypBC8QvnlFw526u2UHvOaB9puRtw+wmQU6OYyxY2Ty+kv2sCqmcOP24i7Q/47g0EOklqY4IKxd/CtpzUDD9yUXtuOOOM9ddd501VDN1h01i5MiR5r777iv0Pky5YJXgwIsFj3D79u2tt5jKEq7KBBdOJcVFhy+++ML6HYHfNegn0Vg4nmg9ykwHg0CmEsvSdSksMOjlYkjUlyoPxZU0A5La2rZta66//nob2RTJwzHBwJ4p5zC0TyZCiB3GzQIGMbFGxLdEcN3F1sV1ndkHZ03i3EZAius1gSg/6ADKpLqcIBFuyoUgOS4pEeyMzX7cqO67774r0wpgU6BeL6MEvF6c6BDAXOQAUYwtgeYchNMRwpRR6927d6H3efrpp+2og8zveIwePdr06dPHe56IMaXeRHSgOYIjUeJk0KGmsfO9MrCLsghmOh4LE8c4F8mSBggc+1wkSZK9+eabbeKMSA22I9sPjy0ChDwMrCNhqA0c1LbnInE0GBHsElo5fqnG9M477xQZ6DJtjvglCKCqONGiYcOGVgQHtXNczud+sDAUB35eujiVxODBg+2SCAzalFMS0YRyWS55hmhXWLt9cRFwIIJjfe9hh+lRZnoQvlOnTi2xtBkJsGeffba9QPI3wMBAArj0MHvnjhWiwUEVweRs+I/ppk2b5nqVRAowG0fTKpLbEMN4g139XwcRX47tss4oi/B0jmsmESxE5ptjhDWawGCtZs2aNkrK6Jmpw7CUsSoOvgsiduzYscWWJQS8ZG3atLGzP61atbL3mR1yUCdWpKdzHL+La3UbNEiQdpDLEdZjOl/Bk0mSKudmfzkxfsfmzZvbhPugd/IUZYeZpyB3jst5JFiIskKkwd9NimLtYYaRMiKYiwZWpLD61vF/0VAHry/ToCXBdCld3Fj8F0c8w+735QKqbmHpE8FB9QWTNOWaJjAI8rfjFeGBREaOfY5hRDH3yfNhoC/yg7oB7xyXlAi+/fbbU3ru1ltvLdtaCZECri6w8w2GvRsYIpgqF84SESYR7JLc6PxIkptrMJMIfquOHTtaywMZ4fHqbxM9IB/A2UXSXSIr36BMETMObFO2LYOtoJUGo2mMq/NJGbeglsoSJU+FU7IUawvRX/J+RH5RLuDJcUmJ4AEDBhSZinL34zWTkAgW2YILOG2x3T6ZSneisPiCwwAjfBrdIH6dj7M4iOYifKnrW9KF0d8gg+lVUXbwYCKCGTxiT/H79oKAGwSCK2spwkkYalGLzCd8I3yxRgSte2uFdNcIlm9LZLuhAhnuQEmeKEyzEfUiIYyKKST8kUlNyaGgwXpRkhDhi++vpHMECVi0MEb8Jnth5D2pM+oiCvvtt19a1j3fQQS7wQURmiCJYMS5qxZCAiQXUCFEeGmXoGpXaEQwbYyFCCJRKIuWyBKBCEYE0m4yKMlLriEJwveVV14psa0pg2KaliB86eiWqlWFEoyulii2kCAOBsKIv11y0HzBeIHdgEoJcUKInIvgZEqUCZFt8DO6GtYkUvltBGGH6Bf1NJ0lItciGLsDCW7U9U1GNFHGjE5uLP5ErFRxUWBQVYj0wUwDiUokoP3www8mKDDIcglxDKCUECeEyCSqDiFCS1TKosWDmqhM/yMKcuULJhrHNh4+fLh5/fXXS/RykbBGE5quXbvaLPCy/h58nhPBZBeri1T6QAC7IvYkLQXFcoMNwm9vCmoNYyFENJAIFqGEVpzz5s2ztzHcRy1hitrATFljB0CkZLOzF1n5RHwRv19//XWxr0WoU/II4YvvK51Zv9hAXIUAWq6GvepH0CBC7zo5Ed0PQhUSdYgTQmQTiWARSvClMpXrfINRaCgRC/YO15acaDDR1UzCZyB88fs68ZkIBBPC99RTT81YwXtVhcgs/i5diOFci2AGtl988YW9XbVqVRsJFkKITCIRLEIH4nfu3Ln2NlPu1JeNIiTHvfnmm/Y2ZccyIYLZlnwGvn9/kmE8aIF6+umn2yQ3Spxl0n7CetFqFRjgEAkWmRPBQUiOwwvsLDfs6/FqRgshRDqRCBahY/78+V4TBkptUfg/qtPVzhec7naTJLoR8X322WdL7OJDRPrCCy+0SW7ZalSBN5Si6u43xhMs0gvlxyjHR9Sf5Dg84Lny1buqIw7VBhZCZAOJYBHqhLgolUWLBeGH1YBawStWrDB//vlnmWwfiMopU6aYsWPH2soTxdX1RXy3bdvWXHDBBbbEWbbFkfN7g6pCZAZ+U6LBWBDYN+jmVLt2bZMLPvnkE68rILaMqA5shRDBQiJYhIovv/zSK+lEmSeqKEQZosGIYAQrEdtdd901pb/n7z766CMrfCdNmmR9lyVFB7t162a6d+9uu/vkgs2bN3veUBLt1CwhczgR7FoV50IEs49Onz7du9+qVausr4MQIj8plQh+/vnnzYgRI2yN1j/++KNIdIGsbiHSDS1eJ06c6N2nKkGUyqIlEsFumhhLRLIimEYbtDFG/CZzPBJtJepLibNcV2FAjBH1BuojU85LZAasJswOAD57ZlayfUzh/V65cqW9zf6d6kBPCCFKS8pXl/vvv9/84x//iPtcLj1lIvrMmDHDTtkCF8rmzZubqOOPxpaUvERnNSJqlDebOXNm0nV9zznnnEBtS1WFyB5169a1xxKJl5TiY8BEQma24Jrhkj+hdevWWftsIYRIWQQ//PDDxXoJhcgEiF83ZYpf9ZRTTsmLAReeYLzBRMFjRTDHIQlkCN5p06bZGqtbt25Nqo3xmWeeadq3b5/Wur7pgJklIsFOpCsqmHmOPPJIK4Jh9uzZWRXBCxcutLMWbtZD1hchRKBFMFOyXEiHDh1qp08pmyREJkHsYYNwdYERcZmqTRs0EPx09sJ6ROIQC5HSqVOnWuFbUmUHfxtjhC8VHni/oLJgwQLvd6YMm8pkZccSwYDD1enNVmMWjuu33nqrUBQ4Hwa2QogQi2BG6ni4evToIQEssgL7GwlxwMX5hBNOyKstTztbOrcheLE6OL9sSTBA7dy5sznrrLNMixYtQiEwVBUi+5QvX94cdthh1paAhYYZBbr/ZRqsF/4kV7XFFkIEXgTfcMMN5rzzzjNPP/206dOnT2bWSoj/guB75ZVXvO2BhzWK3eH80TG6xNE4wC1uAJAMbBuah2AX6dSpkxXQYWHjxo3W3uEGO7SNFtmBurxEZV293jZt2lhxnEn8XmAqQoRhkCaEyHMRzDQsZZT69u1rRo4caVtb+gvZcyJ76qmn0r2eIo+T4WiMwX5F1zCmbqMEHtjPPvuskOhdvXp1Su+BNQTRQoSc7P4wCd/YiL9L5qNahURR9mDQQSSW3wBbBLYU7CiZAnuP8yHXqlUro58lhBBpE8GURnMXJwqcs8QiESzSwc8//2yjUkSkGGhFIRmOJCCOmQ8++MAKXrrfkfSWCvhkqeaA4CW5jTJiYd8usVUh1CAjNwlyrlU1CXKZFKZ+LzBRYHm/hRC5oFQFOIurDhGFi7HIPexjEyZMsJFBRDAXyrB1kXJJbPhcWYj40vmtNFC7l3JWJLXdfffdNpOedre0vY3CMUe03x8ZDHLyXlShPTbNMiiVxm/Bvso+l25IrnYVQJhVPPjgg9P+GUIIkRER7O/sI0SmIEq6ZMkSe7tmzZq2MUbQRRwil8UJXy72pQUxQjLbIYccYg499FArSCiF5ryzUYPt5gbXskLkBgZThx9+uHn55Ze95hknn3xyRqPAxx9/fMa9x0IIkTYRHHQxIsLP77//biZPnlwoGc7vO881rq2vswOxuChmaXB+Z0QvgpeFdrb+CK+bpgbqBfP6KKGqEMGAfY8Oclh0aLeN3SadiahYnNy+XLVqVft5QgiRK8rUj5Rps02bNhV5nAu4EKXltdde86KdJMLtscceOduYRCdJ4kHoEuH9+OOPbdIQQri0UJMVTy8CAOHLbR4rDuwPyXaOCxvUpXVRf5L88qUGdBCheQr7I2XSSNpkvyc6nC78M4nHHntsoAa3Qoj8o1Qi+K677jL33ntv3Cx2oleu2L0QqYIY4gIMRKBOPPHErG7ENWvWFIrwsiDSSgteXhLXSDI66KCD7P90QUs1EWinnXayC62RsVlEqWsjVgiHqgTkHkrsuWNwzpw5toZwOnzndH10yY9UMOF9hRAiVCJ42LBhpn///vakGKULscg9JMGNHz/e268o+cWUaSY/jyYUTPuy4EP+5ptvSv1+RLUoM+UXvESxK1Qo04SLB8liiGAidAiKoLU8Li1UyHBIBOce9jNqNDMgpZoJ/9NxMB1RYHdsEwWOcr1vIUQ4SPnqTJMM4AKMdxMxTIYvmfDUmsxGu00RTYg6LVu2zOsgRckm9rF0Qf1TIruIXUQv1gYS2kpLs2bNrNhl+pj/qZmdyQs7lgi8yEA0mO6NUbNC1KlTJ9erJIyxFgj3u1AurawimBkWjjk3O0K0WQghQieCSWpA+NLtB5ECRKVuu+028/jjj1s/pxCpghh9/fXXvfunnnpqmWqH4lVfuHChjTKyzyJ+KctU2tkLynYhdt1CBQPsCdnE7wuOigiWFSKYEJGnUyPefI4hBpAl+daLg8omrhHKUUcdFZlZDCFEnolgapMCQsD5xDi5XXfddeaOO+4w11xzjYSwSBkuuEzzA8liRJ6S9Za7DldO8LJgc9i6dWuZfLxO8FLHlCniXNfjbdSokXe7LOXXgoSsEMEEaw/HIR0bOY5o7tK6detSD3D9Pn+avAghRChFMNEAfIlE6UhuYLqaExxF+2HWrFmZWE8RYfDh+hNmOnXqlPC12G4QuU7w8v93331Xps9v0KCBrdTgavLi661UqZIJGmwbItLMvGAbKa3IDwqyQgTfEkEEl9kTagZT07c0szPvvPOON6AlGa5y5coZWFshhMiCCMazhwimMgQdhhAhRAjcyVEnOJEKzCJMnDjRu48ARuzFA3sDyXJlAXG73377eYKX//EfhwWiwYhgRMXKlSuzbsnIlBWCyLsIFnRopB41xx0DFv6nZGEqYKdAQAMJoiTECSFEaEUwF6uvvvrKRu7oJsSF7M8//7TRAqaLO3bsmJk1FZGEQZRrJUxGenHF8xl0MU1LIf9kYH8keY19FuGLrYHb2B3CCiLYRc2JBofZFywrRPAhgQ3xC4jZVEUwM4NcH6Bly5YZrfYihBAZF8EDBw405513ntlzzz3tqJ4M+1dffdUKDmq63nPPPSmvhMhPGDhNmzbNu8/+U5zvligukSm/eHIQZWKfdIKXBVtD1GYm/MlxP/30k4mKFaJu3bq5XiURB0r8ERHGhvTll1/aWQgsOckmp7777rv2NjOF6jYqhAgaKRu8EBZEgJs2bWqnrekzT/kbkh9odcsJMxUeeeQRm4lMhICFyIO/wgSJeFdeeaWNgJFRTBkq/iZeeS1sGYgeyrThX/N3s2Mdu3fv7jUd4HZZmiCIskO5L9qouigwkd6ScJFcEtbOP/98c/fdd9v9hdmJN954wzZxufDCC21ST9QEsPMvO+tRmEWwrBDhgH3N3zHOWRuKA6sO5dAee+wxL9kV2xGlNIUQIkhUKIuXc968eTYy0LZt21KvAOL2zjvvtNPWMGLECCuyKWnF1Fvfvn1tkfVRo0bZigEIncsvv9yKAV7nBDA97mni8eCDD9qIIevmT+I499xzbUb9lClT7P1evXpZIYyIF7mPAjOASab6wj//+U8zZMiQtDWgCBtk1xM1xQpB23LaN4dxW8gKER4YUHLeRdxSJYIZm3jtjskVee+996xQdlWEoHz58qZVq1ZZXmshhCiZUl09ifgiIukm5NokUzOYqB61glNJXurcuXOh+4MGDbKRXuc/Q+D26NHDRnaBzyXCQMMDJ4IRyn369DE33nij9z5+rySeNsQv7+ladT7xxBM26swUH9PoIrsQuXVlvhjQYHNIhjAngqULBo6IYAaiHIPJRNCDhKwQ4YIZFepiE91ldg1POsLYDWZ/+OEHa3sguu9qATs4tkl2TdZCIYQQgbZD4AGmkQEXX06ArvkAkbzvv//etr0tLZR8GjNmjM0odh2FqCk5adIkO/XLZxEVRkARjQAy5Ik+ULUCIU6UDO+Zv1QbQhrx5O9VzxQfj9ENSeQ2CtymTZuc1+ANqy/4xx9/NGGOAqsqRDiItUQ4ywMzb//5z3+sMHYCmBk4LG6XXXaZufrqq0OdvCmEiDYpR4IHDx5ss/Nr165tp2MdXbt2tc+5RIhUL4qIXvxj1BtGSOM9hqFDh5pLLrnERr+Y9uUE++STT3oF17/99lv7/4ABA2xSHu1rR44caYUVdWQ5AROhjteOlcecJzUeZDW7zGZwLXa5ACTbyCFZ8e+WfIDfzNX25TcgChxve+bbdkkWSrq5AShRuHTui9kAweQGzxznOpayQ1mOJyK6JDAS/GCfu/322wvlXAA5IlSAcAEG95lBR+cZbRvtM9E7npK9rlQoTeFzonYkI/nLWTlLQWmSdfhbLoxMk44bN87aHyjSzgUSEUzkgWgwyVNvv/229QQjBLBduOjDpZdeahOigKQpIo1PP/209Y9CvEijK+uWCP6Wahix4HcrS0vfWNhZaDri/HNRB2uKK3NGdJ7IfzzybbskCwNF9nsOcgYTfv9l0ME3unjxYm8AxFR7Otdf+0zmtg2WCAQw+MsUIo45jqnI4rzCYdontc9o22ifid7xlOw5KGUR7CoqMN3lhwQdf7Q0FUhkc4lxCGuSLx544AFz//33m5tuuslGhl0XMT4XwUzUFxHsGh24yLGDKhLuhM1J2tWi9UMku7jSTCTa9evXz7vPd2MqGhHiOuSlAzdi4j2jLvb4TZjC52JJtjhRo0QDinzaLqnCgJCIOi2j2TZUTgkDHLtOKFExIJ3HEWifydy24VjFWuY6hiJ6saDtvPPOobYzaZ/RttE+E73jKTY/IW0imBJo+HCdDcHhun5hkygrRGixIRBtYIkVSWxY9wWpGMFUHQlufvANd+jQwd7GasGJm/bOTNcBPmIe4yReXCY+SyzYMtKdkc93Ygljpn8qEMl3F0wGMSW1J86X7ZIqDMaIArMtsfSExXdJWTz3+2NdysTvqn0mM9uGv2EWjoEsJTKjlKSqfUbbRvtMtI6nZD+7QmmiAVgTunXr5j3GiXHYsGH24nbUUUel9H5EehGrXNSJapEYN2PGDDtlTt1gktyuu+46G+ki+oVNAs/vfffdZ/+ez+T52267zU7XcWGlzNqiRYvMiy++6EWFKaGGt5jKEq7KxEknnaTKEFkEqwwiCKjlTAc3UTrwyDuIrIdBBDOLRPIsqEFGOCEIkmoteCGECCopi2DKkSGCqRLhIjoIS6K3RGyvueaalN4PmwL1ekm4ILKA3QEB7GoPI4qxJSC66VqEEKaMWu/evb334DNJqmPdeA1ieOrUqYVKR40ePdqWUWvXrp2936VLF/PQQw+l+vVFGaCyh4PBjaK76RHBzvYTdFQVQgghRKhFMK2SEY94Zf2VE7AN4OF1pc2S5amnnir2eSJGRJlLghrB/jrBsRC9oOGGyA0MdpwIwivkbCmidGA74pjDFhSWznFqkCGEECJIlMqwQf1HGlVQIQJxQ3IZlga8uULEA4uLK4tFFDhexymRPMzCcLzRcARvO0uQPZqyQgghhAgapXYtcwG+6KKL0rs2IpJgUaENtqsl6m9aIkqPE8HA/0EWwbJCCCGECBopF7t99NFHzWmnnWYbVvihDTGP87wQsV5gV82DJifbbbedNlAaaNiwYWg6x8kKIYQQIvQimCQ4yqHtu+++hR4nGW3ChAnm8ccfT+f6iZDDNPiHH35obyN+U60eIhLjtx8FWQTLCiGEECISIti1u0X0+qFwOsTWDxb5DSXtXPFsajKHpalDGMD+4JpNYIdwnuugISuEEEKISIhgSpHB6tWrCz3u7vsrRoj8hrrPNCgBmmJghRDpTY5zpdI2bdpU5JgMCrJCCCGEiIQIdlOwd9xxR6HIE/f9zwvxzjvv2I5/QDJculvkisL1gl2SXJCQFUIIIURkRHCrVq2s+KW+75577mlOOeUU+z/3iUzxvBC///67mTNnjt0QNMWgvrTIv6YZskIIIYSIjAimRTFF+mHx4sXm5ZdfNt98840VxjzO80LMmjXLs8a0aNEi0OW7wkzQI8GyQgghhIiMCN5rr73MuHHjbMcqhK9b6tSpYx8nKizyG3zj7777rr1NK22aY4jMULlyZdsNEegc50rRBc0KQUMdFiGEECLUzTI6duxolixZYoXOzz//bFsbU/rKRYhFfoMNgkQtOPjggz2RJjJD48aNbUMS/Nccj0Hx5SsKLIQQIpId4xC8rVu3Tu/aiNDDrMDcuXPtbXnEs2eJmDdvnmeJkAgWQggh0myHWLZsmbn22mttTeAqVarYhdv4gHlOCOpEr1mzxm6IPfbYw9pmRPZ8wUFpmrFu3TpZIYQQQkQjEkzTA9oi4/MDVx5t4cKFdqE6BJ3kjjnmmMytrQg8H3/8sXf7kEMOyem65JMIJurOMRkUEfzZZ595tw844ICcrosQQghR6kjwihUrzOmnn24FsEuEc7j7PIdIXrlyZTJvKSIInlTnA6VF8j777JPrVcoLsCa5pDM8wa42cy759NNPvdsSwUIIIUIrgh988EGbeIPYRdjce++9ZsKECWb8+PH29r777mtfx2t4rchPFixY4HUUxCZDlziRXUsE1SFybU3i811EGn+yqkIIIYQIrQh+7bXX7HTrBRdcYD7//HPTt29f06VLF3PyySfb20T/eA6RPHny5MyvtQgkskLkjiD5gj/44APvdsuWLXO6LkIIIUSZRDBNMfytkePxr3/9y/5P4wyRf2zYsMF89dVX9nb16tXNrrvumutVyrsyaUFomvHXX3+ZTz75xN6uWLGiad68ec7WRQghhCizCHZT3MVl+teqVcv+77qEifwC4eMaNSB8mDkQ2aN+/fq2PXWuI8FYYmiZDfvvv7/Zfvvtc7YuQgghRJlFsBO406ZNS/ga95x7rchfKwQNMkR2QQAjhGHVqlVes5Js8/7773u3aZcthBBChFoE4+vD79uzZ0+bEBcLj1100UU2+icPYP5BIpRLxmJanhbaIreWCFooZxsSY7/++mt7u2bNmrLECCGECH+d4AsvvNAKXcqfUSoNz2eTJk3sc99//71tjoBIRgQjlEV+4TygoNrAuaNhw4bebdqaN2vWLKuf/+GHHxaKAssSI4QQIvSR4M6dO5uuXbt69YGJ+CB8WFzpNDjzzDPNSSedlNk1FoECH7ATweXLlzcHHnhgrlcpb9ltt9282x999FGhet7Z2A+cCEb8ajAkhBAiMm2TR48ebcuhkfEN/qYZ1IPt16+feeaZZzK3piKQUA1k/fr19vaee+5pKleunOtVyltq1KjhCeFffvnFq+qSDbBBuG6Se++9t9lpp52y9tlCCCFERtsmk3hDY4z+/fvbJLgffvjBiuBddtnFtGnTRglxeYo/Ie7QQw/N6boIYw4//HBP/L733ntZs0QoIU4IIURkRbCD6g9nnXVWZtZGhApK57k2yZTC2muvvXK9SnkP3Rt33HFH89tvv9nfhvrNVapUyeh24bO++OILe5vP0n4ghBAiUnYIIWKhe+CWLVvs7YMOOsirUytyB7+BK02GT9ffvS2TswFbt261t/EC4w0XQgghgo5EsCg1qg0cTA477DCvMgOWCNfEJBNgifILbdUGFkIIERYkgkWpoCye855ikdl55521JQOUILfHHnt4v5NrZ50JyA1YsWKFvU2r7OK6SgohhBBBQiJYlArKornqIHSIU03Y4CXIOebOnZuxz1FCnBBCiLAiESxSBvHrb5ChNsnBw1+mbOHChV75snQnRn722Wf29nbbbWf233//tH+GEEIIEQgRvGnTJtO6dWtbEm3RokUZWykRbGjJ658CZ/pdBIty5cp5/lwGLf6IbbpAAP/5559eYiT1woUQQohIimDKYHExnTFjhtc2WeQfdCNzKAocXFq2bOnZVEheS3eCnD8hjs8SQgghIm2HcO1Qly1blon1EQHnr7/+Mp9++qlXjktT4MGlWrVq1hYB69at82r5poOff/7ZLFmyxN6uX7++adiwYdreWwghhAikCL7rrrus/+/SSy+VEM5DvvzyS7Nx40avMQOzAyK4HHHEEd5tyqVlKgqsxEghhBCRF8HnnHOO9Ru+9dZbpnHjxqZevXrWF+qW3XbbLaX3e+SRR8wBBxxgqlatahcu2q+99lqhblRXXnmladSokRVcRLb4Gz/HH3+8vQj7l7PPPrvQaygV1b17d5ssxMLtTCQL5VNtYDcrIIILpdKqV6/uDWB+/fXXtMwGuP2A2YDmzZuX+T2FEEKIbJNyiy+mQF3Uh4SbVatWec9xP9WIEOL2zjvvNM2aNbP3R4wYYU4++WRbfYBIY9++fc306dPNqFGjrA/5jTfeMJdffrlp0KCBfZ3jkksuMbfffrt3PzZCee6555qlS5eaKVOm2Pu9evWyQvjll19OdRPkLSRGUmkAaM3ratGK4MLxSLk0N7CkXFrHjh3L9J7YKtxswH777Wd22GGHtKyrEEIIEfgSaYhdVyPW3Xb3U6Vz5872ooygYhk0aJAVWK626Zw5c0yPHj1stBcRjHg98MADzYcffljofbgQE5V2iysPBQg3xO+TTz5pI80sTzzxhHnllVdsdEwkx7x582wU0FUDYEZABJ9DDz3U+604btxvWFqUECeEECIvI8GZbMG6detW88ILL9gok/MyHn300WbSpEmmZ8+eNvpLZQo6YD3wwAOF/nb06NE2Wly3bl3ToUMHc9ttt5kqVap4QhpRTDtZB9ExHps9e7bZc889464P5Z9cCShYv369/R8RUVYhEfu93RJkEFBusIMITuc2CPN2yQWpbBtmRZhVoaTZhg0b7P9YkEoDtiIGjuwH2Cx22WWXjO8HqaB9RttG+4yOJ51ndA7+K8nrUsoiOBPMnz/fil6K7xMFHj9+vNlnn33sc0OHDrVWB2wT+A+JaBHRRRw7unXrZpo2bWojwJ9//rnp37+/jVpOnTrVy2SvU6dOkc/lMZ5LxJAhQ8zAgQOLPI5POZ1RUC7cv//+u71dvnx5E0Twkn7zzTfedmOAwXbIJGHYLrki1W2DbcGVtnv77betf780vPvuu2bz5s3eezpbRFDQPqNto31Gx5POMzoH/5akPqlQWoWNnQCLwerVq82sWbNsFJYoMVHY2rVrp/R+RGIpu0Wi2rhx46z9YebMmVYII4KxRhANJurEBRxPMGWZTjjhBPv3iGQHF+bdd9/dTgGTvOPq2MbzKpfkYUZM9+vXr1AkmGRAhDpLunDRPN4zqGKP36BixYr2NhF1F2XPJGHYLrki1W3DccEx88svv9hmJwjoeAPD4uD4XrBggd0POG6OOuqotB4H6UD7jLaN9hkdTzrP6Bz8d5KuhZRF8JYtW0y7du2sGPWLSOwIJK3dd9995uqrr07pPek05RLjEK94DrE73H///eamm26ykeFOnTrZ55nGRTDfc889ngiOBeHLhfrrr7+2t4kQuw5nfkjqwz6RiG233dYusRCRZkknCBmWdL9vOuB3Zpu7yhv8RtlazyBvl1yT6rZhtgUfvLO2dOnSJaXPw4ZEvWH2gb322svUrFnTBBHtM9o22md0POk8k9/n4GQ/O+U5fUQuUdrYRLiLLrrIPuYvb1ZaeB+8uAhulljrARu3OJVPtIq/I/LlLv5cvP2tY6mZymNHHnlkmdc36nz//fdeaS0GK/6kQxEe/IMXrBEcI6mghDghhBBRImURTMSXSFCfPn0KPX7MMcd40aJUINL7zjvvWKGFN/jmm2+2yW/4fKkbfNxxx5nrrrvOPvbdd9+Z4cOHm5EjR5pTTz3V/v3ixYttaTQiW7zH5MmTTdeuXW3tUqZrgdrC7du3t7YJpvVZuH3SSSclTIoT/49qA0cDKqhQWcWVu8M3nywMGPHbAxYIIsFCCCFEXolglxzlr8kLbmq0uESzeGBToF4vYrRNmzY2QovXuG3btvb5MWPGmBYtWlhRjEeYmsKUUevdu7dnpZg2bZo58cQT7XsgzrFrvPnmm4W8koh3WvzyHAu2imeeeSbVr593kATlxBLWEKoMiPDir5CSTAc5ZgAmTpxo/ud//sfz29IkRfYUIYQQYSdlw4azJsQm4zhxnOrF8amnnir2efy8w4YNS/g8iWrYM0qiRo0aNnlPpAbWEqp2AIOIeB5pER5ILsUmtHz5ctv4ZtmyZbb0YCwkz3FcUU7Nbz2qXLlyoVbMQgghRN6IYNoiMy2KJcEfzXXJcOoiFi38TUnwlIpwg5WJaPCECRO8aLCzFuHFJ5kU65Eb1DqYcWFGBtsTA0ohhBAi70TwaaedZr27V155pVcZomHDhl6lCJ4X0YCSdU4MIXyoxSzCDxVTSGAl+RS/N355uioS+SVC7IfIL956EkjVHlkIIURei+Brr73WjB071l40nQh206XUIu3bt2/61zIPIXMfK8KiRYvMmWeemZMWxQgkVwUEH2hxNZVFeNhuu+1sxz+iwAhhmsI4y4vf409SKr+7qw8thBBC5LUIJjJEc4xbb73VTqlihaDWLlOqdFejRasoOww0XEIaQoQGINkE8es6jLl1ENEBS4RLjPMLYDz2xx9/vE2AzMXASwghhMgWpapkXL16dfPggw/aRWQGqlc4EUxENtsi+IcffrDNRIAWu/KBRgvakOPvp8SgKyN47LHH2t9aEX8hhBD5QKnbeVAPeOrUqWblypW2/SolzZQUlz6ow8q0NVE6EhFPOeWUrFZm8EeBlRAXTc4//3xra6I6BFVYhBBCiHyiQmmmya+66irz6KOPFuoaR/SI2r1EhxVJKjv4MGls4Hyb+INJaMqWH5k2yYDwpjSaiB5Yl7K1TwkhhBBBI2XT37/+9S/z8MMP22Q4RLBbuP/II4+YwYMHZ2ZN8xC/QPF3bcs0qg0shBBCiKiTsgh+8sknvQjhOeecY2688Ub7P/cRw48//ngm1jMvadKkiefFpX4rrWuzgRLihBBCCBF1UrZD4AHG7vDSSy+ZDh06eI9Td7RTp072eZEe2M5Eg2kBzQADiwJlqzIJQhu/t0uAJFFKCCGEEMLkeyTYVSmIFWNklrsscxFeS8Qnn3yi2sBCCCGEiDwpi+D+/ftbkfTCCy8UevzFF1+0/99yyy3pWzthatWqZXbZZRe7JejmtWzZsoxtFX5XtUkWQgghRD6Qsh3ijTfeMLVr1zY9e/Y0TzzxhBVoS5YsMXPmzLHtk19++WW7uOn8p556KhPrnVc0b97cbmMXDaakVSb48ccfPTuLagMLIYQQIsqkLIJHjBhhxS1RQ4QvC3CfKCXP+5EILju0uGVgsXXrVusL7tixY0a6efmjwOoQJ4QQQogoUyol5eoD+0ukxd731xAWZWOHHXawzTNg/fr15ptvvslIbWDXoa5SpUqqDSyEEEKISJNyJHj69OmZWRNRLERmqd/rLBHp7s73xRdfmE2bNtnbNMegW50QQgghRFRJWQRnukSXiA+RYDp8IVRpo0wXuXS2UVZtYCGEEELkEymLYMfUqVNtkhyJVHXq1DHt2rUzbdu2Te/aCY8KFSrYNspz5841mzdvtkI4Xb5dagN/+eWXXm3g3XbbTVteCCGEEJEmZRH8xx9/mK5du5rJkycXevy+++6zCVuUTtNUeuZqBiOCXeQ2XSKYZDvn4eY9SXwUQgghhIgypaoT/OqrrxZJgmNBGN98882ZWVNhy9HVrFnTbonFixebtWvXpr02sKpCCCGEECIfSFkEjxkzxkYKaeKAIH700Uft/9xHUD377LOZWVPhtVEG10a5rCxdutSsWLHCqw3sRLYQQgghRJSpUBr/KFC39rDDDvMe79y5sznyyCO950VmQATjx3ZVIkhULIt9QQlxQgghhMhHUo4EUz4L9t1330KPu/suUikyA5HaJk2a2Ns///xzmdoo//XXX+aTTz6xtytWrKjawEIIIYTIG1IWwf/6179st7IHHnig0ONDhw61TRYGDx6czvUTcfAPNIgGl5aFCxeqNrAQQggh8pKU7RCI3J122snceuut1g+88847mx9//NH89NNPpnbt2mbAgAHea5mmnzZtWrrXOe854IADzKRJk2wkF19wp06dStVG2Z8Qd+ihh+b9dhVCCCFE/pCyCJ45c6bnQWUq3j8dv2rVKvu8S9xSqa3MtlGmVvCGDRvMV1995bVVThb+ztUGrlatmmoDCyGEECKvSD18+F+BW9Iigm2JwAv8999/29uqDSyEEEKIfCPlSPB3332XmTURKUHkl4jw77//bhYsWGCbmCTbpES1gYUQQgiR71QoTcMGEYw2ygcddJCZPXu22bJli5k/f75p0aJFiX/3yy+/mFdeecVWlgAqTVDjWQghhBAinyiVHUIEg+bNmydtiSBSTEc/2lt/8cUX3uPHHHNMRtdRCCGEECK0keDy5csnfI4peKLDJ598srnxxhtt5QiRHajMQRSX6O63335r1qxZY6pXr17E+oBAfu2118z69eu9x/mdqCrh6j4LIYQQQuQTSUWCXaJbvAS4TZs22SoDd999tzniiCNs1QGRmzbKrvGFg9J1Dz/8sHn++ec9AYyNonXr1ubaa6+1dgohhBBCiHwkaTtEcRUfnCBGDN97773pWjeRYpUIRDC/AwORsWPHmgcffNAsWbKkUFc/xG/79u3Ntttuq+0rhBBCiLylQlkrQmzcuNE2bOjfv7+NPE6cOLFQwwyRWWrUqGGaNm1qf6MVK1bYJhofffSR9QA76tata7p06WJ23313/RxCCCGEEMmK4JIqQuyzzz625uz5559vFi9erA2bg2iwG6i8++67hfzabdu2NUceeWSxvm4hhBBCiHwjbdUh6tevb///888/U/q7Rx55xLYBrlq1ql3wFZPE5fjtt9/MlVdeaRo1amS23357s/fee9u/iQdWgA4dOliv7IQJEwo9R9JY9+7dbUIYC7fXrl1rogDbD6+vg+9/2GGHmeuvv95Wf5AAFkIIIYQoY53gRLz44ov2/5o1a6b0d4jbO++80zRr1szeHzFihK00gb8VD2vfvn3N9OnTzahRo2xN2zfeeMNcfvnlpkGDBvZ1fu6///6ErZrPPfdcs3TpUjNlyhR7v1evXlYIv/zyyybsMDg49thjzVtvvWV23XVX07lzZ9OwYcNcr5YQQgghRLhFcM+ePRM+R8eyzz77zCbFuQhkKiDY/AwaNMhGeufOnWtF8Jw5c0yPHj3M8ccf74nXxx57zHz44YeFRPC8efNsDdwPPvjAi0o7Fi5caMUv7+nW74knnrBRZ9Z7zz33NGGHZLdWrVop4U0IIYQQIl0iePjw4QkjrH4rAq/p06ePKS1bt241L7zwgk22Q6DC0UcfbZO9EOJEf2fMmGG++uor88ADDxQS4uecc4556KGHTL169Yq8L0IaC4RfoB9++OH2MTquJRLBWDv89g5XZuyvv/6yS7rge7ulLGB7SOd65Zp0bZcoom2j7aJ9RseSzjM6/waNrQG5bierhSqko0Sam5KnVjDRyFSh5S+il4oGO+64oxk/frxNtoOhQ4eaSy65xNom8L2WK1fOPPnkk1YcO7BMkPwVa49w0CK4Tp06RR7nMdc+OB5DhgwxAwcOLPI4PmXWI12wsyDkQf5dbRftMzqWMoHOM9ou2md0LOXLeea3335Lnwi+7bbbSuwYRxWCVP3ADiKxlFkjUW3cuHHW/jBz5kwrhBHB2BiIBvM5b7/9tvUEY3k44YQT7ON4YWMbRcQSL5LtoteJoOxbv379CkWCGzdubIU6S7pwIybeUyJY20X7jI6lTKDzjLaL9hkdS/lynvn777+zI4LTQaVKlbzEuEMPPdT6erE7kOh200032cgwLX5dJQQE8z333GNFMAKYsmzVqlUr9J6nn366rYyAfQKLBDV0Y1m1apWtoZsIGkrEaypBRNpfjSEdsLOwpPt9w462i7aN9hkdTzrP6BwcNHRtCva2SfazA6m4iNDixd2yZYtdYq0HbFyn8m+88UZz8cUXF3p+//33N//+97+9pDusFuvWrTPvv/++admypX3svffes49hoxBCCCGEEPlFzkUwkV5q+2IzoN3vmDFjbPSWag7UDT7uuOPMddddZz3H2CGwSYwcOdJWggCivPGS4XbeeWfbSQ2oLUz1BLzFVJZwVSZOOumkSFSGEEIIIYQQIRPB2BSo17t8+XJbrQG7AwIYjzEgivHmduvWzfz6669WCFNGrXfv3il9zujRo23linbt2tn7tBGmmoQQQgghhMg/ci6Cn3rqqWKfJ8o7bNiwlN4zXiWLGjVq2IYbQgghhBBCpK/OlxBCCCGEECFBIlgIIYQQQuQdEsFCCCGEECLvyLknOEw4r7Frn5zO9n50N6Hsm+oEa7ton9GxlAl0ntF20T6jYylfzjPr/6vTSup2LBGcApRwA8q5CSGEEEKIYOs2Ko8lYpuCkmSy8GBks2zZMlOlSpVi2y2nimvH/OOPP9rayELbRfuMjqV0o/OMtov2GR1L+XKeKSgosAK4QYMGRRqu+VEkOAXYkI0aNTKZgh1GIljbRfuMjqVMovOMtov2GR1L+XCe2amYCLBDiXFCCCGEECLvkAgWQgghhBB5h0RwANh2223NbbfdZv8X2i7aZ3Qs6Tyj828Q0LVJ2yXq+4wS44QQQgghRN6hSLAQQgghhMg7JIKFEEIIIUTeIREsQgc1midMmJDr1RAhQvuM0D4jdJ4RsUgEp4EhQ4aYFi1a2CYaderUMaeccor58ssvixRuHjBggC3cvP3225vjjz/eLFiwwHv+119/NVdddZXZc889zQ477GB23nln06dPH7Nu3bpC77NmzRrTvXt3W/+Ohdtr1641YeaCCy6w2yyfyOY+M2jQIHPkkUfa11SrVs1EAe0zmdtnvv/+e3PRRReZpk2b2vfYbbfdbKLL5s2bTZjRPpPZ80yXLl3sc9ttt52pX7++vTbRXCrMaJ+pk9F9xvHnn3+agw46yAYrPv30U5NNJILTwMyZM80VV1xh5s6da6ZOnWp7Z7dr185s3LjRe83dd99t7rvvPvPQQw+ZDz74wNSrV8+0bdvWa8XMyYLlnnvuMfPnzzfDhw83U6ZMsRcjP+eee67dSXiOhducbES4yOY+g3jp2rWrueyyy7L+PUX49plFixbZ7piPPfaYvbD9+9//No8++qi56aab9HOGjGyeZ1q1amXGjh1rBdO4cePM4sWLzRlnnJH17yzCs884rr/+eiuocwJtk0V6WblyJa2oC2bOnGnv//333wX16tUruPPOO73X/PHHHwU77bRTwaOPPprwfcaOHVtQqVKlgi1bttj7X3zxhX3fuXPneq+ZM2eOfWzRokWh/Rl79OhRcPLJJ9vbu+yyS8G///3vQs8feOCBBbfddpt3n+87fvz4giiRqX3Gz7Bhw+zfRwHtM9nZZxx33313QdOmTQvCjPaZ7O4zEydOLNhmm20KNm/eXBBWtM8UZHyfmTx5csFee+1VsGDBAvs5n3zySUE2USQ4A7iQf40aNez/3333nfn555/taMpBDb3jjjvOzJ49u9j3oe1ghQr/1916zpw51gJx2GGHea85/PDD7WPFvY/I331GRJds7jO8xn2OCC/Z2meYDh89erS1YVWsWDHt30NEY59ZsWKFueSSS8wzzzxjbRO5QCI4zRCo7Nevnzn66KPNfvvtZx9jh4G6desWei333XOxrF692txxxx3m0ksv9R7jtfhHY+GxRO8j8nufEdEkm/sM09oPPvig6d27d1q/g4jePnPDDTeYypUrm5o1a5offvjBTJw4MSPfRYR/nykoKLCea84rhx56qMkVEsFp5sorrzSfffaZee6554o8h+nbDztB7GOwfv1606lTJ7PPPvvYhJTi3qO49xHhINP7jIge2dpn8PW1b9/eesovvvjiNH4DEcV95rrrrjOffPKJeeONN0z58uXN+eefb99LhJMrM7jPMLDmuf79+5tcIhGcRsiGnDRpkpk+fbpp1KiR9zimcYgdJa1cubLIaApjORedHXfc0YwfP77QVBLvw/RBLKtWrSryPmGlXLlyRU6aW7ZsMVEl0/tMPqB9JjP7DAKYZKcjjjjCPP744yZKaJ/JzD5Tq1Yts8cee9gkqTFjxpjJkyfbBKsooH0mvfvMW2+9ZfcNrBRYJJo1a2YfJyrco0cPky0kgtMAoo0R00svvWR/WEoL+eE+Ow6Zlv6MfbIw8Uw5GBXhs6lUqZIVRpSa8cPFCF/N+++/7z323nvv2cf87xNmateubZYvX15om+BBihrZ2mfyAe0z6d9nfvrpJ1v26OCDDzbDhg2zAiBKaJ/J/HnGBTMofxUFtM+kd58ZOnSomTdvnq1wxcKACZ5//nlb1jNbKHsmDVBO5Nlnn7X+J+q+uhESCWvU0GOK4JprrjGDBw82u+++u124jRGckmduxMQO8/vvv5tRo0bZHYjFHXxMLe299952VIWRnPJF0KtXL3PSSSfZenxRoHXr1racSufOnU316tXNLbfcYr971MjWPgN480hU4f+tW7d6dRgZeTNCDzvaZ9K7zxABRgBT25MSR8w0xUaBwo72mfTuMwRmWPCOct7+9ttvza233mprTBO8iQLaZ3ZK6z7D+cWPuxaxz/hnRTNOVmtRRBQ2Y7yFklQOyopQ5ovSIttuu23BscceWzB//nzv+enTpyd8n++++8573erVqwu6detWUKVKFbtwe82aNQVhpnv37gWnn366vb1u3bqCM888s6Bq1aoFjRs3Lhg+fHgkS6Rlc5+hzE+81/D3YUX7TOb2Gd4v0WvCjPaZzO0zn332WUGrVq0KatSoYd+jSZMmBb179y5YunRpQZjRPmMyem3yw+O5KJG2Df9kT3ILURSi20QlKbwtRDJonxGpon1GaJ8RsUTL6CVCBS2gX331VTNjxgxzwgkn5Hp1RAjQPiO0zwidZ0S6kCdY5IyePXvalov/+Mc/zMknn6xfQmifETrPiJyja1P+IDuEEEIIIYTIO2SHEEIIIYQQeYdEsBBCCCGEyDskgoUQQgghRN4hESyEEEIIIfIOiWAhhBBCCJF3SAQLIUQIaNKkiW1ZmsxC63F3mxbIQgghiiIRLIQQQggh8g7VCRZCiJBGhpcsWWJvT58+XRFfIYRIEUWChRAiYnz//fdx7RADBgzwHn/sscfM9ddfb+rWrWuqVKliTj/9dPPzzz+bZcuWmbPOOstUrVrVNGzY0PTu3dv89ttvhd6/oKDAWi6OPfZYU61aNVOpUiUryq+44gr7HkIIEQbUNlkIIfKQW265xaxatcq7/9JLL9nI8rp168w333xjH9uwYYMVy4he/gdun3vuuWbMmDGF3o+/ffjhh8348ePN7NmzrSgWQoggo0iwEELkIYjZ999/3yxfvtzssssu9rGPPvrIbNy40Xz22Wfm66+/NpUrV7aPjxw50r7eiWUngC+44AL793/88Yd59tln7WPcv+6663L2vYQQIlkkgoUQIg+56KKLTIsWLUy9evXMYYcd5j1+4YUXmv333980a9bM/g+I3BUrVtjbRHodWCLq169vtttuOxsddkyZMiWr30UIIUqDRLAQQuQhiFzH9ttv791u2rSpd3vbbbf1biOEwYnh4sBD/Oeff6ZxbYUQIv1IBAshRB5SoUKFlB53kEjneO6556xNInb5+++/CwloIYQIIhLBQgghkubUU0/1bvfv39/MnDnTRolJqJsxY4bp2bOnrRIhhBBBR9UhhBBCJM1pp51m/b8kwlGKLV5Huh49emiLCiECj0SwEEKIpKHG8KhRo0yHDh3M008/bT799FNbSq1WrVq2ysQJJ5xgzj77bG1RIUTgUcc4IYQQQgiRd8gTLIQQQggh8g6JYCGEEEIIkXdIBAshhBBCiLxDIlgIIYQQQuQdEsFCCCGEECLvkAgWQgghhBB5h0SwEEIIIYTIOySChRBCCCFE3iERLIQQQggh8g6JYCGEEEIIkXdIBAshhBBCiLxDIlgIIYQQQph8438BloJOmifJO04AAAAASUVORK5CYII=", "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# -- Plotting upper ocean heat content time series for Pacific Ocean -- #\n", "plt.figure(figsize=(8, 4))\n", "plt.grid(True, lw=1, alpha=0.2)\n", "\n", "(1E-21 * H_upper).plot(lw=2, color='0.1', alpha=0.6)\n", "(1E-21 * H_upper).rolling(time_counter=12, center=True).mean().plot(lw=3, color='0.1')\n", "\n", "# Axes labels:\n", "plt.title('eORCA1 JRA55v1: Pacific Ocean\\nUpper 500m Ocean Heat Content', fontdict={'size':12, 'weight':'bold'})\n", "plt.xlabel('Time', fontdict={'size':12, 'weight':'bold'})\n", "plt.ylabel('Upper Ocean Heat Content (ZJ)', fontdict={'size':12, 'weight':'bold'})" ] } ], "metadata": { "kernelspec": { "display_name": "dev", "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.13.12" } }, "nbformat": 4, "nbformat_minor": 2 }