tssic_cmems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://www.eurocean.org/np4/file/157/1.png\" alt=\"drawing\" width=\"200\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sea ice concentration service\n",
    "## Introduction\n",
    "The aim of this service is to generate a monthly sea ice concentration field for a user defined time period.\n",
    "\n",
    "The service uses CMEMS product “Svalbard ice chart” provided by the Meteorological Institute of Norway. The product has a resolution of 1 km and is provided on a daily basis. A detailed Product User Manual is available at https://marine.copernicus.eu/documents/PUM/CMEMS-SI-PUM-011-002.pdf. \n",
    "\n",
    "This notebook has been created as a part of INTAROS project. This project has received funding from the European Union Horizon 2020 research and innovation programme under grant agreement No 727890.\n",
    "\n",
    "Torill Hamre, NERSC. 2020-05-25."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "This service requires the following libraries installed in your notebook environment:\n",
    "\n",
    "* netcdf4, https://github.com/Unidata/netcdf4-python\n",
    "* matplotlib, https://matplotlib.org/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the libraries needed\n",
    "from netCDF4 import Dataset\n",
    "from datetime import datetime\n",
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing for output\n",
    "\n",
    "The notebook will generate output in NetCDF format, using the same map projection and coordinates as the CMEMS \"Svalbard ice chart\" product. Hence, we start with obtaining information about an instance of this CMEMS product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CMEMS product variables: {'time': <class 'netCDF4._netCDF4.Variable'>\n",
      "int32 time(time)\n",
      "    long_name: reference time of sea ice file\n",
      "    units: seconds since 1981-01-01 00:00:00\n",
      "    calendar: standard\n",
      "unlimited dimensions: \n",
      "current shape = (1,)\n",
      "filling off, 'yc': <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 yc(yc)\n",
      "    axis: Y\n",
      "    long_name: y-coordinate in Cartesian system\n",
      "    units: m\n",
      "unlimited dimensions: \n",
      "current shape = (2980,)\n",
      "filling off, 'xc': <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 xc(xc)\n",
      "    axis: X\n",
      "    long_name: x-coordinate in Cartesian system\n",
      "    units: m\n",
      "unlimited dimensions: \n",
      "current shape = (3812,)\n",
      "filling off, 'crs': <class 'netCDF4._netCDF4.Variable'>\n",
      "|S1 crs(maxStrlen64)\n",
      "    grid_mapping_name: polar_stereographic\n",
      "    straight_vertical_longitude_from_pole: 0.0\n",
      "    latitude_of_projection_origin: 90.0\n",
      "    standard_parallel: 90.0\n",
      "    false_easting: 0.0\n",
      "    false_northing: 0.0\n",
      "    proj4_string: +proj=stere lon_0=0.0 lat_ts=90.0 lat_0=90.0 a=6371000.0 b=6371000.0\n",
      "unlimited dimensions: \n",
      "current shape = (64,)\n",
      "filling off, 'lat': <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 lat(yc, xc)\n",
      "    long_name: latitude\n",
      "    units: degrees_north\n",
      "unlimited dimensions: \n",
      "current shape = (2980, 3812)\n",
      "filling off, 'lon': <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 lon(yc, xc)\n",
      "    long_name: longitude\n",
      "    units: degrees_east\n",
      "unlimited dimensions: \n",
      "current shape = (2980, 3812)\n",
      "filling off, 'ice_concentration': <class 'netCDF4._netCDF4.Variable'>\n",
      "int16 ice_concentration(time, yc, xc)\n",
      "    long_name: sea ice concentration\n",
      "    standard_name: sea_ice_area_fraction\n",
      "    units: %\n",
      "    coordinates: lon lat\n",
      "    grid_mapping: crs\n",
      "    source: met.no\n",
      "    _FillValue: -99\n",
      "    scale_factor: 1.0\n",
      "    add_offset: 0.0\n",
      "unlimited dimensions: \n",
      "current shape = (1, 2980, 3812)\n",
      "filling off, 'concentration_range': <class 'netCDF4._netCDF4.Variable'>\n",
      "int16 concentration_range(time, yc, xc)\n",
      "    long_name: concentration range\n",
      "    units: %\n",
      "    coordinates: lon lat\n",
      "    grid_mapping: crs\n",
      "    _FillValue: -99\n",
      "    comments: Range of the analyized ice consentration value\n",
      "unlimited dimensions: \n",
      "current shape = (1, 2980, 3812)\n",
      "filling off}\n",
      "Variable ice concentration shape: (1, 2980, 3812)\n",
      "Variable ice concentration has num lines, columns: 2980 3812\n"
     ]
    }
   ],
   "source": [
    "# Get dimensions and positions of CMEMS SIC data for Svalbard region\n",
    "filename = 'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806011500.nc'\n",
    "sic_dataset = Dataset(filename, 'r')\n",
    "print('CMEMS product variables:', sic_dataset.variables)\n",
    "print('Variable ice concentration shape:', sic_dataset.variables['ice_concentration'].shape)\n",
    "yc = sic_dataset.variables['ice_concentration'].shape[1]\n",
    "xc = sic_dataset.variables['ice_concentration'].shape[2]\n",
    "print('Variable ice concentration has num lines, columns:', yc, xc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate sea ice concentration statistics for given time period\n",
    "\n",
    "After selcting a month and year, the notebook obtains all CMEMS Svalbard sea ice concentration (SIC) products for this period and calculates the mean SIC values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing day 1 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806011500.nc\n",
      "Processing day 2 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806021500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806021500.nc'\n",
      "No SIC data for day 2 of month 06 of year 2018\n",
      "Processing day 3 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806031500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806031500.nc'\n",
      "No SIC data for day 3 of month 06 of year 2018\n",
      "Processing day 4 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806041500.nc\n",
      "Processing day 5 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806051500.nc\n",
      "Processing day 6 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806061500.nc\n",
      "Processing day 7 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806071500.nc\n",
      "Processing day 8 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806081500.nc\n",
      "Processing day 9 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806091500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806091500.nc'\n",
      "No SIC data for day 9 of month 06 of year 2018\n",
      "Processing day 10 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806101500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806101500.nc'\n",
      "No SIC data for day 10 of month 06 of year 2018\n",
      "Processing day 11 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806111500.nc\n",
      "Processing day 12 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806121500.nc\n",
      "Processing day 13 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806131500.nc\n",
      "Processing day 14 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806141500.nc\n",
      "Processing day 15 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806151500.nc\n",
      "Processing day 16 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806161500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806161500.nc'\n",
      "No SIC data for day 16 of month 06 of year 2018\n",
      "Processing day 17 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806171500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806171500.nc'\n",
      "No SIC data for day 17 of month 06 of year 2018\n",
      "Processing day 18 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806181500.nc\n",
      "Processing day 19 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806191500.nc\n",
      "Processing day 20 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806201500.nc\n",
      "Processing day 21 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806211500.nc\n",
      "Processing day 22 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806221500.nc\n",
      "Processing day 23 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806231500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806231500.nc'\n",
      "No SIC data for day 23 of month 06 of year 2018\n",
      "Processing day 24 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806241500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806241500.nc'\n",
      "No SIC data for day 24 of month 06 of year 2018\n",
      "Processing day 25 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806251500.nc\n",
      "Processing day 26 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806261500.nc\n",
      "Processing day 27 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806271500.nc\n",
      "Processing day 28 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806281500.nc\n",
      "Processing day 29 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806291500.nc\n",
      "Processing day 30 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806301500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806301500.nc'\n",
      "No SIC data for day 30 of month 06 of year 2018\n",
      "Processing day 31 of the month\n",
      "File name https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806311500.nc\n",
      "OS error: [Errno -90] NetCDF: file not found: b'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/2018/06/ice_conc_svalbard_201806311500.nc'\n",
      "No SIC data for day 31 of month 06 of year 2018\n",
      "\n",
      "Number of days with data 21 of month 06 of year 2018\n",
      "SIC statistics data added up: [[    0     0     0 ... -2079 -2079 -2079]\n",
      " [    0     0     0 ... -2079 -2079 -2079]\n",
      " [    0     0     0 ... -2079 -2079 -2079]\n",
      " ...\n",
      " [-2079 -2079 -2079 ... -2079 -2079 -2079]\n",
      " [-2079 -2079 -2079 ... -2079 -2079 -2079]\n",
      " [-2079 -2079 -2079 ... -2079 -2079 -2079]]\n",
      "SIC statistics data sic_stats: [[  0   0   0 ... -99 -99 -99]\n",
      " [  0   0   0 ... -99 -99 -99]\n",
      " [  0   0   0 ... -99 -99 -99]\n",
      " ...\n",
      " [-99 -99 -99 ... -99 -99 -99]\n",
      " [-99 -99 -99 ... -99 -99 -99]\n",
      " [-99 -99 -99 ... -99 -99 -99]]\n"
     ]
    }
   ],
   "source": [
    "# Loop thru all data for chosen year, month\n",
    "stats_year = '2018'\n",
    "stats_month = '06'\n",
    "num_data_days = 0\n",
    "sic_stats = np.zeros((yc, xc), dtype=np.int16)\n",
    "sic_data = np.zeros((yc, xc), dtype=np.int16)\n",
    "for stats_day in range(31):\n",
    "  print('Processing day', stats_day+1, 'of the month')\n",
    "  str_day = str.format(\"{:02d}\", stats_day+1)\n",
    "  filename = 'https://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-metno-svalbard/{}/{}/ice_conc_svalbard_{}{}{}1500.nc'.format(stats_year,stats_month,stats_year,stats_month,str_day)\n",
    "  print('File name', filename)\n",
    "  try:\n",
    "    file = Dataset(filename, 'r')\n",
    "    num_data_days = num_data_days + 1\n",
    "    timestep = 0\n",
    "    sic_data = np.int16(file.variables['ice_concentration'][timestep, :, :].data)\n",
    "    sic_stats = np.add(sic_stats, sic_data) \n",
    "    file.close()\n",
    "  except OSError as err:\n",
    "    print('OS error: {0}'.format(err))\n",
    "    print('No SIC data for day', stats_day+1, 'of month', stats_month, 'of year', stats_year)\n",
    "  except IOError as err:\n",
    "    print('IO error: {0}'.format(err))\n",
    "    print('No SIC data for day', stats_day+1, 'of month', stats_month, 'of year', stats_year)\n",
    "\n",
    "print('\\nNumber of days with data', str(num_data_days), 'of month', stats_month, 'of year', stats_year)\n",
    "print('SIC statistics data added up:', str(sic_stats))\n",
    "\n",
    "# Calculate statistics, whole %s for SIC\n",
    "with np.nditer(sic_stats, op_flags=['readwrite']) as it:\n",
    "  for x in it:\n",
    "    x[...] = x / num_data_days\n",
    "print('SIC statistics data sic_stats:', str(sic_stats))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store generated monthly mean sea ice concentration\n",
    "\n",
    "The calculated monthly means for SIC is stored in NetCDF, using the same map projection and coordinates as the input CMEMS product.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lat_var VALUES: [[59.92085  59.923893 59.926933 ... 54.375866 54.370758 54.365646]\n",
      " [59.928665 59.93171  59.934753 ... 54.382217 54.37711  54.371998]\n",
      " [59.936478 59.939526 59.94257  ... 54.38857  54.383457 54.378345]\n",
      " ...\n",
      " [78.695206 78.70398  78.71276  ... 67.15177  67.14316  67.134544]\n",
      " [78.6967   78.70548  78.714264 ... 67.15248  67.143875 67.13526 ]\n",
      " [78.6982   78.70698  78.71576  ... 67.15319  67.144585 67.13597 ]]\n",
      "lat_var SHAPE: (2980, 3812)\n",
      "lon_var VALUES: [[-21.288599 -21.273005 -21.257406 ...  38.812492  38.823395  38.834293]\n",
      " [-21.294678 -21.27908  -21.263477 ...  38.821266  38.83217   38.84307 ]\n",
      " [-21.300758 -21.285156 -21.26955  ...  38.83004   38.840946  38.85185 ]\n",
      " ...\n",
      " [-80.275154 -80.26748  -80.25978  ...  85.25425   85.25609   85.25793 ]\n",
      " [-80.31994  -80.312294 -80.304634 ...  85.27642   85.27825   85.28008 ]\n",
      " [-80.36474  -80.357124 -80.3495   ...  85.2986    85.30042   85.302246]]\n",
      "lon_var SHAPE: (2980, 3812)\n",
      "Mean sea ice concentration data is stored in file tssic-mean.nc\n"
     ]
    }
   ],
   "source": [
    "# Store statistics in NetCDF file\n",
    "outfile = 'tssic-mean.nc'\n",
    "dset = Dataset(outfile,'w', format='NETCDF3_CLASSIC')\n",
    "dset.createDimension('yc',size=yc)\n",
    "dset.createDimension('xc',size=xc)\n",
    "month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']\n",
    "month_index = int(stats_month) - 1\n",
    "dset.setncattr('title', 'Monthly mean sea ice concentration ' + month_names[month_index] + ' ' + stats_year)\n",
    "license_stmt = 'This dataset is made freely available by NERSC. \\\n",
    "User must display the following citation in any publication or product using the data: \\\n",
    "\"These data were produced by NERSC as part of the INTAROS project. \\\n",
    "This project has received funding from the European Union Horizon 2020 research and innovation programme \\\n",
    "under grant agreement No 727890.\"'\n",
    "dset.setncattr('license', license_stmt)\n",
    "dset.setncattr('summary', 'Mean monthly sea ice concentration from CMEMS daily products for the Svalbard region.')\n",
    "dset.setncattr('spatial_resolution_latitude',  1000.0)\n",
    "dset.setncattr('spatial_resolution_longitude',  1000.0)\n",
    "dset.setncattr('southernmost_latitude',  54.365646)\n",
    "dset.setncattr('northernmost_latitude',  88.10237)\n",
    "dset.setncattr('westernmost_longitude',  -80.36474)\n",
    "dset.setncattr('asternmost_longitude',  85.302246)\n",
    "dset.setncattr('time_coverage_start', stats_year + '-' + stats_month + '-01')\n",
    "\n",
    "crs_var = dset.createVariable('crs', 'c') #char\n",
    "crs_var.setncattr('grid_mapping_name', 'polar_stereographic')\n",
    "crs_var.setncattr('straight_vertical_longitude_from_pole', 0.0)\n",
    "crs_var.setncattr('latitude_of_projection_origin', 90.0)\n",
    "crs_var.setncattr('standard_parallel', 90.0)\n",
    "crs_var.setncattr('false_easting', 0.0)\n",
    "crs_var.setncattr('false_northing', 0.0)\n",
    "crs_var.setncattr('proj4_string', '+proj=stere lon_0=0.0 lat_ts=90.0 lat_0=90.0 a=6371000.0 b=6371000.0')\n",
    "\n",
    "lat_var = dset.createVariable('lat', np.float32, dimensions=('yc', 'xc'))\n",
    "lat_var[:] = sic_dataset.variables['lat'][:,:].data\n",
    "lat_var.setncattr('long_name', 'latitude')\n",
    "lat_var.setncattr('units', 'degrees_north')\n",
    "print('lat_var VALUES:', str(lat_var[:]))\n",
    "print('lat_var SHAPE:', str(lat_var.shape))\n",
    "\n",
    "lon_var = dset.createVariable('lon', np.float32, dimensions=('yc', 'xc'))\n",
    "lon_var[:] = sic_dataset.variables['lon'][:,:].data\n",
    "lon_var.setncattr('long_name', 'longitude')\n",
    "lon_var.setncattr('units', 'degrees_east')\n",
    "print('lon_var VALUES:', str(lon_var[:]))\n",
    "print('lon_var SHAPE:', str(lon_var.shape))\n",
    "\n",
    "yc_var = dset.createVariable('yc', np.float32, dimensions=('yc'))\n",
    "yc_var[:]= sic_dataset.variables['yc'][:].data\n",
    "yc_var.setncattr('axis', 'Y')\n",
    "yc_var.setncattr('long_name', 'y-coordinate in Cartesian system')\n",
    "yc_var.setncattr('units', 'm')\n",
    "yc_var.setncattr('_CoordinateAxisType', 'GeoY')\n",
    "\n",
    "xc_var = dset.createVariable('xc', np.float32, dimensions=('xc'))\n",
    "xc_var[:]= xc_pos = sic_dataset.variables['xc'][:].data\n",
    "xc_var.setncattr('axis', 'X')\n",
    "xc_var.setncattr('long_name', 'x-coordinate in Cartesian system')\n",
    "xc_var.setncattr('units', 'm')\n",
    "xc_var.setncattr('_CoordinateAxisType', 'GeoX')\n",
    "\n",
    "sic_var = dset.createVariable('sic', np.int16, dimensions=('yc', 'xc'), fill_value=-99)\n",
    "sic_var.setncattr('long_name', 'sea ice concentration')\n",
    "sic_var.setncattr('standard_name', 'sea_ice_area_fraction')\n",
    "sic_var.setncattr('units', '%')\n",
    "sic_var.setncattr('coordinates', 'lon lat')\n",
    "sic_var.setncattr('grid_mapping', 'crs')\n",
    "sic_var.setncattr('source', 'INTAROS')\n",
    "sic_var.setncattr('scale_factor', np.float32(1.0))\n",
    "sic_var.setncattr('add_offset', np.float32(0.0))\n",
    "sic_var[:] = np.int16(sic_stats)\n",
    "dset.close()\n",
    "print('Mean sea ice concentration data is stored in file', outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate a plot of the mean SIC field\n",
    "\n",
    "Using matplotlib, the notebook generates a PNG image of the computed mean sea ice concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Graphics of mean sea ice concentration is stored in file tssic-mean.png\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAE/CAYAAAD2T5JXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9e5gU1Z3//zoGYWC4jAMLYyM4Iq0DwTgC35WL+4WA+1sHFYyIm7C4wUs2EMxlo1kTNTpec1l1o5HFZL2wATVBcYWwjLuREfcHKFlAFJYhzgiNyDCwiANy6XHJnu8fn3O6TldX9/T03Jl6P089XXWq6tSp6qp61+eutNaECBEiRIgQIZqGM9p7ACFChAgRIkRnREigIUKECBEiRA4ICTREiBAhQoTIASGBhggRIkSIEDkgJNAQIUKECBEiB4QEGiJEiBAhQuSAkEDbGUoprZQanmbdXKXUurYeU2eHUqpCKfXV9h7H6YbwuoYIkYyQQA2UUjGl1GdKqQG+9q2G5Ipb4BhrlVK3NLefEJmhtS7TWv9ze4+jvdAS95lSqlwptdRta63rGnSstoBS6qtKqc1KqaNKqY+UUj9VSnVz1hcqpf5FKXVcKbVHKTXbWXe2UmqlUqo26P1g9v2NUuqQmZ5XSvVtu7ML0RYICTQZu4Gv2AWl1EVAz/YbTogQLQ+XJLo4egHfAQYAlwJTgdud9QuBz4BBwF8Bi5RSnzfr/hd4DZiZpu8HgbOAYcD5po/ylh1+iHaH1jqcJBtTDLgb+E+n7RHgLkADxaatH/Ar4L+BPWafM8y6ucA6s98nCCGXmXUPAX8E4sAx4EnTroF5QLXZZyGg3P7M/ELgUd+Yfwt8J835aOAbpt9PgQeQB/kt4CiwDOjubH8VsBWoBzYAX3DWfR/4wPSzA/iSsy7tOacZ1x3APtPXH4Cppv0M5zgfm/EVOvu9BNQBR4D/AD6f4RhrgVuc5a8BVc74R5v2CLDc/Je7gW9l6LMn8Kj5z4+Yc+5p1k0H/stcu7XACN99dTvwntnvN0Ces36Gue5Hzblf4dxnzwD7zfV6EPhcM++zBeZ+2G3aHgf2mmNvBv7MtF+BEMf/mD7e9V9X83/dba7HQeSZ6GfWFZvjfRX4EDgE3JXh2pYDS337dgv6PzOde2PXLYt3wHeB35r5fHMNLnDWLwF+7NunG877wWmvAL7hLC8A/q2933Ph1LJTuw+go0zmRXc58lIfAXzOvFzOJZlAfwWsAPqYh/194Gazbq556XzN7D8fqMUjxMSLwDmuBlYBBcBQ5GV+hdOfJdA/NX1Zsh4AnAAGpTkfDawE+gKfBxqANcgXcT+ESL5qth1tXoKXmnF/1VyPHmb9LIRszgD+EjgOnJ3NOfvGdKG5phGzXAycb+a/A7wNnAP0AH4BvOjse5O55j2AnwFbM/yXietsxr4P+D+AAoab//QMhDTuAbqb67IL+Is0fS40/Q425znBjOUCcz3+HDgT+DugBvNxYq7j7831K0SIfJ7znx4x+55h+i4x61411yAfGGj6+Hoz77PfmTFY4p8D9EdI4DbkAyXPrCvHkFqa63qTOc9hQG/gFWCJ879q4J+QD4+LkftvRJprmzgW2RFopnNPe92yeAe8iiFI4BLgpG/97RiCddrSEehVwGpECj0LqCTNx244dd6p3QfQUSY8Ar0b+BHyFf479wExD2wDMNLZ7+vAWjM/F6hx1vUy+xaZ5XQvtsuc5WXA953+1jnrqoA/N/O3AqsznI8GJjrLm4E7nOVHgZ+Z+UXAA779/wBMStP3VmBGNufs2284QtSXA2f61lVhpFGzfLZ5UXYL6KfAHKNfmvG5L9x/A74dsM2lwIe+th8AzwVsewZwErg4YN0PgWW+bfcBk537ao6z/qfAU2b+F8A/BPQ5yNxnPZ22rwBvNPM+m9LIM/CJPUcaJ9A1JEtYF9r/C48Ez3HW/x74cprjJo5FdgQaeO6NXbdGzv1G4CNggFn+M6DOt83XMM+605aOQCPA64iq93+Rd0n3xsYRTp1rCm2gqVgCzEYe1F/51g1ApJU9TtseRHKwqLMzWusTZrZ3I8esc+ZPZNj+nxGpAfO7pJF+DzjzJwOW7XHOBW5TStXbCRiCvARQSv21caay60Yh1yJl/JnOWWtdg0ia5cBBpdSvlVIRZwz/4hyjClFFDlJKfU4p9WOl1AdKqaMIKeEbQzoMQVSjfpwLRHznfCfyEvZjAJCXpp8Izv2gtf5fRMoOvCdI/n8zje1MYL8ztl8gElVKn024z/a6C0qp25RSVUqpI+YY/cjumoLvvM18N5KvX7b3dVOR7tyzuW4pUEpdA/wYUQUfMs3HEO2Ni76IGSAbvIRop/qY/T4A2txRKkTrIiRQH7TWexC7yjRELeXiEPKVfa7TNhSROLLqvpnDWwrMUEpdjKiZX21mfxZ7gYe01gXO1Etr/aJS6lxEFXcr0F9rXQBsR9ShTYbW+gWt9WV4qvGfOGMo840hT2u9D/mgmYFIrv0QKYUsx7AXsf0Gte/2Ha+P1npawLaHEJtiUD+1OPeDUkohxJjNPZFpbA2INGTH1ldr/fmAbYOQ7j5LtCul/gyxR18PnGX+1yN417SxezXpvJHn4BTJH2m54Lj57eW0FWW5b5Ovm1LqCuT+vlprvc1Z9T7QTSkVddouRmzd2eBi4Bda6+Na62PAU8g7JcRphJBAg3Ezou467jZqrf+IqFgfUkr1MeTyXbL/sjyA2Ixygtb6I+A/Eclzudb6ZK59+fBPwDyl1KVKkK+UulIp1QexJWnENotS6kZEAm0ylFIXKqWmKKV6IIR0EpEyQV4wD5lrilLqT5RSM8y6PsiL8WPkxfpwEw77NHC7UmqMObfh5hi/B44qpe5QSvU0Uu4opdT/8XdgpMpngceUUhGz7XhzHsuAK5VSU5VSZyK2xAbEEasxPAPcaPY9Qyk1WClVorXeD/w78KhSqq9Zd75SalKW55zNfdYHIbz/RojiHpIlrgNAsVIq3TviReBvlVLnKaV6I//Jb7TWp7IcYyC01v+NfHzMMdf5JoI/MoL2bdJ1U0pNAZ4HZmqtf+/r6zjyAX2/eR4mIh9xS5z98xA7OEAPs2zxn8At5t7qCfwN8G425xGi8yAk0ABorT/QWm9Ks/qbyFfyLsQb8AXk5ZoNHgeuU0p9opR6Isfh/TNwEY2rb7OGOdevAU8idrAaRIWN1noHYi99C3mpXgSsz/FQPRBV2SFEDTcQUZuCXJuVwL8rpT5FHIouNet+hagI9yHOT2834dxeQjxTX0DUb68i3r1/BK4GShGNwyGEbPul6ep2YBvyYjyMSM5naK3/gKjTf276uBqRZj7LYmy/R2xv/4BIf2/iSXV/jZgLdiD/ycuIXTgbZHOf/RviKfo+cm3jJKt4XzK/HyultgTs/yxyD/4Hcv3iyLORK1yJ92vA95APps+T3ceIRVOu2w+R/3u1UuqYmSqc9d9AnKAOIh8M87XWrgR6ElH1Auw0yxY3IZqSj5D7dhjmmQpx+sB6roXoJFBK/V9E4i02klGIEJ0aSqnHkI+R77T3WEKEaApCCbQTwagIvw08HZJniNMBSqkC4C+AdBqfECE6LEIC7SRQSo1AAvXPRuIgQ4To1FBKXYV4p25EbMkhQgRCKfWsUuqgUmq701aolPqdUqra/J5l2pVS6gmlVI1S6j2l1OhWG1eowg0RIkSIEB0ZxnR1DPiV1nqUafspcFhr/WOl1PcRb/I7lFLTEHv8NMSP4nGt9aXp+m4OQgk0RIgQIUJ0aGit/wNx3nMxA3GqxPxe47T/SgveBgqUUtk64DUJIYGGCBEiRIjOiEEmdMmGMNmEGYNJ9ij/iOTEJi2GDlGVYcCAAbq4uLhV+q7dvJnImDGt0rd7DJCAyRAh0iEvw7p4Dvu46PE5M3MmkkvIRif2Qizn/4PcoHVIkMdn0OBEOfewN+9nsr7huNdnwx+9MdrxxIEzx+TR65M4B876E/rwKb12xyUCNp68X48zkTxBhdLBp3n59DkhB6+tgsiYMZzgv+i1O07tYZKe19rNm4lcBCe6e1ei16dxSapZBxSBNqsOvQN/MmYMHNzsnVgd8IXWff47IzZv3nxIa/0nLd3vcKUSqaGagv2SoMJ9DH6ptf5ljsMISrDSKrbKDkGgxcXFbNrUOk545UpR3kp9u8cAL2gxRIggBFZNd1DjbFeTacOAvoY70avqOwjTTUXyBVUh9WKWQfXtEB0F1Ru97aP+G3co6DVmTH6lmYM/bJpAhP2MvqFKUqXnmcmm2dgL1bUQ/RNgommbA5XTv8CUrW+xs/Rcfq32wObNlH8IrITyW6F80yZ5bo1/RrlSlN8Xh8t8A3iKRBTxltIRjGaHLCz03p+2vxDJUErtaXyrpuMEkhy8qSiHuNZ6bBN3O6CUOltrvd+oaA+a9o+QbGAW5yBPQovjtFfhlodOUiE6ABojz6Zsn2mdmgr6YdD3IGH7T5kVK5O3i17qTYAk4rOT6ScpiZ1/DIVwZf9K+ighz4pamUDI2RJ0NOLsNAqWTp9JLRF2lp5L8fE9/FH/gAH6Zlk/HQbomxMfpBbb9RKuu2cJO6dIjon4pcAaqH4Jqi+RafTOKm7hSaayCnWOTkwtjvvN2AbnlMnytIdCpLKmTjliJVI5CvO7wmn/a+ONOw44YlW9LY0OIYGGCBEiVfK0RBnUlg5JkuNGsBxYfiuUzwrYIYqQZnXy/onjFXrzNYdTl/2oTvedPwq4DMayiU2Mpfj4Hmrzi3jgpR+xatYUDg/JYzE30p3PEoR6nZMh80HupmTvHlgDeb4xRucAD8HT875J6cS3JB8UZJ8WP1vcr+AeQ8pfxyPRfQ5RD1aUm2vQFT/eFWJFaPF+lXoRmAwMUEp9BNyLZDVbppS6Gak7a+/w1YgHbg0iFN/YCkMCOgGBuqqcjo6NhGrcEM1HDUKUwwsNaeXQh923eiM8v1Fe5j9RijvcG9SSp52vNhKsISgVBV2d2qd/rDWGMPzjtGQata+1Kug+pYFenGBH/ghG762CCETYz6t8ibk8R+G8OIefyiPGeUxgA2uZDEDJXXsSJM96YKIhTgeVE8fzLZ7gZsZlfZ2ahHt08vw9CIm60ug+TXk2fQ1WycR7msBKoC0NrfVX0qyaGrCtRgqYtzo6PoF2EvIMESIIuahua+x0GK4ozLx9zWGxaepaUK669CFghmRK7wlwg+IOSzhDne0cggQgAuoGM1+bun64fzyHkwnfjgk8yTl6HPTjMr5h+XUMi1YkjWH0+ipGD62Shb1QuD5OIbI8mipPQo4Ci0nYU+NGPb0jfwQABdRTylbuvvlOlvDX7FlRQs5YmKqiXXSrJINOeift00kq5/Js+z8NyRNaTwLtqOjQBNqZpM8QIYLgOgbltP/hANIyUBHPxphEngB3QfROOPbQA0znBVhYJYQ4D88uGkCeQJK7hZoH+ikywiXPxLidef02qHGIku1RZ8VQ02bnG0MPEk5D5bfC95/zyNOPb/EET8z4VhadZoHHM68O31EeWksC7ajo0Oea5IXXSW7SUI0bIghWSstlvyBJ1BKmrhW1K4AR2gSr5aeYGP/CNTyy4HZe2HuzR56WLK1kZ8nUkqcNgwHUNKg25kg/UQaRe5AHsfbVz3ntMFzhErWduQf09ORt1VRgDsSnwo+dktx5eyFWUpy0bT0FnEV9wAiagIUKHodFvg+MFOkzRApCCTREVuhMpB6i/ZGrBGrxmlGLJhFpvvGU/dBb1tXQ8CHkmYqpt/KMuFOsdvaL4BFl1Cx/iBQPs2RqifVDb7d0krAL12abLhTnNZ8d1ZWe9TzfujuBcfD2xIsZtzO1nGYxMWKmvvo7lKaQ52JdxSClmO9/Vl0V7QLfuseB9zUY1aytEB4+740jlEDbATYRgQs/QQURVnuSmCsdhwiRCbnaQYPwmiUoo9pNhJoYu6UC8haTTJJW6nTJyfWWdaVDK3X5Qliic9JLoemQzXmnqJ596/TDoO6H2okRbim5hXP4ZtKzN10LiVpnoyCc1PMThLnoVmmb757f48nP8KJqEuSZQrwhMiKUQDsQ/OTUEaW+cq1TxnkAGNQ+wwnRBeDaVa2XrFoi84Gxm25bJnum3w75YfKi3+tVr05eDgprsQ5OQXjtMJRlINCEeroWrl1Zwerp0/ij/gGVTAAkZrQXr3OCXnxCQdp+nuBbsAC+u3BRom1RtY9EnXYI1bW5IpRA2wH788agLpJ5vQ3UjRLkEwRLoh2RTC1ihAQawoO1f0YzkIVF2jjKALhqXUukbthJgkytOvYloD+iqrXwL7twnXz841ya6nWbC6pNcaroqMzbMRRK2coENgBQy9lspZQC6vmM7gAp6lvrhXvujJ18d+GihPRp4bdxhmg+Qgm0PRAHdspsNjnzQ7VpiM6EJPL0Sz3OS7wp5OnC9dR1JVBdbVSk/iTN/X3LpUidiw/Jzhs2KhmMdBYEVGEI0tpu/fZPC0ukrnpYuRF+6+HW9c8AsGXBCK5aWMnGBZfyGd3ZbWygLp7A88Dds6KERbeKOnZRgGp2UYAGqaN+nHd0hATaXjhlb1jF8v3TeC8LjuzIUmiIEOCoKNOlxbNOO0PIOVtnDSRsogmVbtRHQEOB46n7Al7dCkusrqRs5/1S6l4ywm/DTUecfqhIgMq3GklEb8h99MIq4nNhM2PZzFjGkJrr1o0BfXTGNwCoyMKuGdo8m4+OQyqtjw6ZC/c9VZH1tkF20qD2pvbTXBxofJMQpzGG47PvVTvTEGdqQTVizWFvAjz1qw1vGeLbYUhAm1/NPIJgFe8Q53g5jzh1XyuJqnnIx4X74fEh7JpVxK4FRWzIH5+2zyf4FufO2JmYblvxj5zU8yl7EuavSJU4i50pRPNgJdCmTp0VHeNjIQ8YLjd11bZift1ELsuGRBuTVF1noKZItemIN0ZoB+3KqCG90JmQ3oYAU5BKJoiaV9uCTr5sPpngersmVMWuhDgCkeDsMffiEedeMwZ3bHadVedayfRj37ID/zizqSiTbt/oRf4G4HrYOeRcutPAQQZRm8L0yfgWT3gLM+TnsQXz+e7KRcxfIdJomdYJqTREyyB0ImoHjPk8eBWHYgSXc2seypXir5DUZuX+Z8+k1QrVwSFaFXciZAaSqMB6xFoCuw3Uo8BxQyIDoealzF0GkmeiAa/EsCXD9ebXEmWQVArQH+LmCyDPVf1uBT72VMVREy7jJ8HmSKWJcV5N4kOCG6BkwR6IwLChdewYMpIxbApU41rynMZqSm7Ykwi/ARH47dhc8iwLn/0WQVezgXZIFW5rEJklT/A5a3wdqbJwf6oUa6dcsbHxTUKcpiibRLIIasnMkFz5ZXjq271IvYiLkFR306H6TaefNMgqvvSgmVxvWjeJgl9MHioZf+Klspi31be+VCY1T8gzV8enTKh409dgSX8hQuyOo5Elzs2klpIsmbyHCoc8bY5hP0LybDm0cTmzdkeHHXtQfGUufSQtI8T4PPBXtRBdB7ilkdYrWAPlvhia5pJomNqva6EsQoIcq2shaklwJbAQdq0totyWJ7xfCYE+h5DTbRomwisPK+446O+5iRhIonIJG0lWz7rE6XjexqNQm18EwLCtdeLBexy0kZxVFaLCjdAsx6fGUPEmlFUD3zXncZCENO3PfxuhlloibHJI9AZ1M9aTIpM0vBEoC2hvivknRNdFhyVQaB6JprvpU9onGmnTkKmfPFsCIYme/iibRKpD0Nd9At5q4DYYhpCn7q9Q38Yjob2yHuAOe59ObsZH5Po07XNJG99pcZBBfFYqCXELqOcEPTnvcJ23wRDgBogaT9+ah3MfZkYcRh4eq86phZFDqqjNl3g3G8JSYGJAZ/M8V62spIJkCT3IRmvRXG/+MBrAQ1dT4XZoAgVS0vm15jHKL0vtvy2OH6Jzo2wSIh35CfQwMB2RLi2THiRRmPm1wzD8XohG8Uh0soK13j2Xos70odFqL/a4fvuoTZLgi/uM5Z9Lvcnq4/525zN2F8Iw6oiXQt5hYCkwCYhD2RyS1KUtgiie81PE+83bCziVyoTge3nkOcPE3mpNtXlmy5x5F8PN5CZoyRbh+yAVoRNRB0ZLqHWb0n/QV2WuhBpKoacnEuQJHllVAz9DwkfWAQ8jJHobsBWefOhmzjf3TtQRUReZMJf5bv8ZPEVTQkBqAxyJqp3xDcFLED+U7JImGESopeD4EeKlouItyK8ntuw8Rh6vokcDqDWiCq322RybhYkBbU+RlNN3P5GE9Amwc/q5wB6ABGFG05Cni78KaMskVYZSZzC6mgTaIZ2IMiHbmzYXonXDWBojz6aMJcTpibIIHjnhzE+BLbNGsOWeEXAjbBkyQrxvDwLVMIENlE2CsvuS+5sfhfm+tmyQyZkoycnHl/zAetnGo95UnyanbC0R6vP78XL+TDYwkX/gbylmNyfy89hdWER8GmxZMiIpk1BzKtCURRCpMwPRb2BCgjxv5xEAdishz6jWRM3zWaGUV6Q8wzGrjac+hM92ruhqTkSdjkCh9W7udMTZUmMJvXJPH5RFCJaQBiJqW4SMtgwxDi+ztKz7Nox+tAquhvht3m7l1RB/J7kNaHKcYhJhRiE6y8zXkprSD49EQSTLGoYnJhCSeppbWM00XufypH1Xc2WCcF/Pn0ItZ9PwYXY5fzOhbBJipwX56ICkax03jlDFxACYzQspfdjrFkSYjRFpY2hMmu3K6GqJFDolgWZCLiTYEscM0XVQFiE4S8I8M/WH0VurRGqjgBeYLetnaSGyoRCfB4/3RupOIvdQXr5mrykY/ROlcnpRJ5GXKx1bic5KoR9CXrXYPC1e5UtJfb3O5WxgAvUUEKOYV7mGegqop4BBHKSeAjYwkV6cZKMxUMTyz4Ulkis3V1S8iXgsp0GeOYdS3km0HVBVfE5VJvdjrp+1c7rzuUrH1Y6WKkQqQgLtJGgPosyEbMcSSqGnASaSTE4gxFnoTZWl44lQywvMZgcjedaSaBWUXy/k2RfggmSStGpH64WbqdxXk+AmTXBItGTvnsQm01jNQQZykIHEKE4UqgZPtbuWyUlTPQW8zuXsJ8L7XMgGJvDDKT+Aaanlz5oCXUtwiEwUWAOXH68kYryZR19ZFbChYLhv/nk89W7QBGIPDTIBhZJnduhKKlylOwAJjR07Vm/alJoQujMine01yAEqdCrqXEhyGLKwxGlQWTqeegoSqeY2MIE3mMxD3M1NRtV4Ow9ygl58iyd4ntk8sPBHUpFoMbAt4HmcrNJ64wZJUtEg9bKbkGEq8D5wide0asgU1vLFhE1xK6XBBzTYwcjEvFWlAjzI3Yy+vwo2QrWpF9pUdak9p+idZsZ6N0eQdII2leBEYCWsWjKFC33Sp4uo845T3/PaZ/79Un6gbmBMwDswE1k+T9tIoP7UolnFpg422+zLPD6l1GatdWr2iWaiVCldmQMj9j9Fq4yntRESaAvC7wqfjddwSKKdA9mS5z/yDUawg/1EqCVCxTPXikcuUus2yXtzpYLpZn6xuU/mpj6PmeygfgINJE+AOXiS51QgX0gzwn5TW/MS6iloMoFOZm2iLifA6q0z4SnQL3lJ7XMmUH8WJiNB20LeahpQDYffzqNwZZzqGal9RecAS+Sans2uRHvd94ZJX9c4pHRZZmJ0vXpbHYPT/+fltY0QaCPkCa1HoJecofSbPZq+X7945yTQTqvC7ShIm8he66Rff5udyrQO1bodHIHkCSnkuZVSaonw4PqHeWbFrVR88CV4NXmXpBffdGd+rpZpofNCVwouUBygCarcIPK0mIqQUL4kZr/8eCUjj1cxgI9TvG8LjLUzCJY8dzCaf+S7bHCy0e8qLUpK4JAreYKXzjCBvT7yjALjoHBhnC3TR9BH9yP6pKzfqGcyT/9WtrlfMZtnE918kbUU/f0uvvL3z0pyCzOV30lGuGreVkUG8oSAXN5Z7tcWUArO/FzTp86Kzqx+7hDINtzFzavrjyFz58PqEB0QQeTpxCL61ba8CvQG7laJQvGUwLPMTqhx06H8VrjvKZFWQbxzy6NAD4mzrMg1dd5e4AI4PD2PgwyihvMTXrk1DKeUrWyllGJixChOqGWtA5GFJc/r1Rhw7tvVXAvALhbAQ6DuwksE30TYYt1qGimhNyqKpDy8IXjfugX96LMALud1Lud1HrtHomrHsokItbzO5USo5Yus5YW9N+c2wNZEEAkGSZTudrYYRibJtI2gFHTrQqzShU61/ZGNSrfMsXeE6t32R9pk7o70OWXrW7xSWsZqplG1cLQ0voyQ5ykN3eQ/t+R5NrvYz7DUDylzb4zfVkm5murdL+9rGKdgBJQdTJ+dKGP4yBA4XJIHQA3nE+M8YpxHMbsTnrUF1CfI06pwYxRTS4QJbJDze2Y0+maSyNOFpClUqUWxmwg1FRhnFvzFu+9CioPnkwgTiqRJymtDXB7hdkCI9ZFxP/S8qJuQTKLNkIUK1t0mlzKMrQUFnNmFWKULnWrbIEX6DCDMhESaKdNJ6HTU7khR3Vrbzl1e06qSKXxGD0/6PAa8BmzXUKSgDrgGiIO6CJ7ZNpv9g19MKy2Ua005JBFU4l65SAVKdQnizKS+vcCbHc4HAAkS9atvLXnWEqE7nyWkUgv1PdB/n+FYqzOsawTRiJE+LcGNI0GgehuolXB4SF5i/JCePC0Kjh9hfO8fMvNSKN9oikrY/m1OYEOk7ZphKEv7pR8dgTgTUEAnVsk2FSGBthHSqnWzIFG/lBISaesirc0TklS3q0qmUM9ZHGQgNQyn6t3Rhjyroci8oYuMGvduYKzmpvtfhH2GJH3wp5NL+fjalpzWz9oLkyq+BKGHJ30m1LcGxexmq+OO+waTOUmvJM9aK4X24kSGgzh4SsNLKrdkBVFfNeC3SRQBV9O8ZkueIAT5496pXdm63Nvs70bnmVoq62faUJuAvMBtig5gv2wRdLFkuKETUQsilyxGjal0g4h3I2E8aZsjikie/SUTTnwIjGQHQCJrT9HFuyQcBaDuKAxApuuAtcBSSSKf9hD34cWLEpB3uTLDvZKO8NPAElA9ZyVUuPUUUPnuVdRTwJYKuMIAACAASURBVFZK2UopMYrZTTENdKdqxeikJO6ZoJobv9qAhK+MM9NAEgRXuD5uRrWblaoqkDxBiHObs+x34tsGLF8qk14N+imp4R2iGehiufxCAm1jpJNEm9pH6MHbynBrcc4j4bRi08jZmplWtRmjmLpnJDQiIX3uNNPLSFL57aA+0PQ8cjjpP7ekueueIj5ULyYN4yP9c6brEUzXIyQGMgAZpU9IK1lZL9tStlJAPT19EubuI3JugzjIeVYiLZYvhGzu2Vyy/STsu5aAq/ESQJh40MMTRZp+QsXJFunU5TMNkb5iVONjRgEXKTa3kzPfaVHhJSTQEK2B1no4Qom0FTAQL/3dseRVtflF7GAkG5iY8FKtJSLxnjFE+uSAOBCd0kKi2zUcAh5UxF8upFxrnmU28eOK2cdfJH5cMWx9XcqLfjYvMLqyitGVVbAwh/PoQWopM0R1C94HQC9OsN9sWNsQ4QS9iG8vTNh2dzSMhLGGsErgvtsb17Y0qxrLVuSDwUbIOFmNCtcHE6cbHpbOyz3dfjO15pXDJJyfxowKDlFLeNNflDo16/nep1neAbxoQzQNIYG2AYJiQoPWN6dvCEm0xWDVobXIC/wST/IEL5TjIAPZwUjPgQggDyRJ3yDTcFBIdJSCBxGHoivi9DxymJsWirSZdwHk5WuYmKqNmPKdt6SE11PJ3rdWuqvBxEseJBjj4PACkdg2MCHRXM9ZAAkHIfkgmAAxOLKpiANHBkJxnALq6cUJRvYQdTWb8kQt3Qiqt+eeb7biTUTarCT5pp6Dl4UIKH/SmW8B4rEkuny7EGm5S6IXKVk2UxDKR8HyZn4kN3f/dod1Imrq1EkRZiJqZbSVV5/74g2djHLDcIw61LUnGvVofK78WukT4OqFa0SqfNlse8j81h2Fbn1lfgDo/aCuBmIwc5sUzJzLcxxkkIS2XKngX809cr2CWcBCL5FAJkkubVWYCF71lRvg7ZKLqWE49RSY0JWzEtJnA915cN89jBi8g7OppXLf5YwYvCOhwr2k31Z2m23rVgxLnNe9l6W/t3V/xWs5xoGmnF8Uzzh5JWwpkdCVegr4D/UWkEZF24zwDktk15pwpfsynEsQoS7fLr8zszy2S5zZ7tMctFYmorE9ld50XuPbpYynqnNmIurE2ufOgbZSybi21LB4d9OQFOtpVbe1eKrPqfJj7Z4D+Fgch15DVLYx4NhRRPI05HnK/O+HFOpsIdFypUDBdD2C0RdVEX8bfphfzOdWA/ZDa5mGyQqqhTht5ZAa5xd3Ph15gtwEez3CGU4Nr3M5Mc6jgPqE9FnFSIoG11K1byQnBveCQ3kwGM7rF5PzJsJ5xIREx8ZFCk3juGOhPtYSVd9MVNTC8FoTtmmyC47eWZU4py/oMt5TFWnJMtfnb6Z5nrY5xJmOpMu3p5LoTLO8XKk2IcQOgy7mhduFTvX0QKaE0rk4JHVVpJCmiyl4qkNHdRs5Xsfr+VP4RxZImr4UVAODRMVZB6K+HYgNzEj8X4sVFduhrBoe4EfJSQnGKSrS6OJrfPNpkzzUIpKzyTxUy9lJKlsgkTChngI+oYCJbGA511HfUMCIi7dINw0RGuLdGdRP9MPnEYPBcHJAL45sKkpz8FZELaLOnpfc/AVdxrXrKySXbQsmFcimD7+ndCCRXpT6TC7fnippzrQfPoMV5bXpx9CusarZoBOrZJuKRm2gSqkhSqk3lFJVSqn/Ukp927SXK6X2KaW2mmmas88PlFI1Sqk/KKX+ojVPoKshyEkiCKE9NBhlkxziGUgqedp1ruON8X6dkr+GWUdeomL9tbBdwSbEwejYZqCv8b7tI6rcIgXdBhGIuZqyFcCHsO6S5PSNy50/zkqZTXLGsY5PQ5CKK8AENlDAJwBJZcpiFFNgAlkAxg/eQEGPemobIhRQT6RHaoKCiWyQmXjbOrxUvwnVLyHORfeIFFpAfSIJxM6J5/J/9fjE9u1FMuXbs9tuprGX2mmmz9HLxgRXm7qw7mTLrXXIj+UuFsaSzdBPAbdprbcopfoAm5VSvzPr/kFr/Yi7sVJqJPBl4PPIo/y6UuoCrfUfW3LgIQQd/mu0gyCJNIPwbWd+q/n9ELmDN8JZoz4mvr7Qs3PuRCTN7ZsRBeNRqDsADDfqW+WpcYNgshnVAHMnIUnkF2i+QHaEWRYhtWA2Zigm8QBA4c44h0vyqOcsQ6LJqfpApNLnmZ0g0iOHCmCwrJvQbwMfOO5A65lAzx4naBh7mHLVv03uPVdtDYjD1Dwkthaxhe5gJO9QyudIbxNtbSQdM0Dq9GPmKM9W6qLafLf4E2skHcv83x3u+e9iKtxGJVCt9X6t9RYz/ylQReLxCsQM4Nda6wat9W7k/v/TlhhsiFS4tQL9X6ShFOqDSzjzzPRtksnTpqFzhK8JD60hvt14k3yEkGcMRJCzOeH6ynw3JR63dRo4yMljZ6D3JydIYJziwHZYbF+e30VIb1zy/5eRSKMEhqckwZFCx7KJrVySsHu6lVZiFDOIg1T9ZjR71nuZEmyiiAa6m7QFxYl18UNnNXLwlkcNjoeuD/UU8MDeH/FlfS7QAeIpt2lvyoCErdTcb9VZ5BB2cx6nrczSnuhCXrhN+lZQShUjZXg3Iu4Ltyql/hpRZt2mtf4EIde3nd0+IoBwlVJ/A/wNwNChHTGjc+dCUFJyP4LCCpoVq9dJEJiarwfQH24p+XmiqYB6mAiTeYOrdlZye8kDgLyc80Yd9kh0HSKJbgchToPJZp2pwKLXFUG+/C+zj78I+aYSy9ua45wNqo65j8DiGTB3nayKXgo12Xz5VCNONf4XaDUwHZFArRS6N07hx3uoL5UMQ6VGxH6VaxI20T8ccZLlxvLoOViSKkxgAyfpxW6KE45E9UcK6HfOgSwG2TqoqIWyKIxeXwUTvdy9zw75CjftfZEv63P5tdrTcaQzP4n6pFMriQbZPTMV9k4gxxy6rYIuJoFmfapKqd7AcuA7WuujSqlFwAOANr+PAjfhS2VpkPLvaq1/CfwSJIyl6UMPEYR2//LuQEirtv020B+mlSxPNEWoTUqqvrRkZiL/a7H1QO2Np8IFX4KFo7C2L3KrH4Rug3hlYpkp8mXiPB3UcxZz9X4WKcVRgIka3pbwj0bPyxcPmZCW70SqlHxIkgTNSli6YGaCOK0d1FZfScJO6Denjg1HJnBhv/fZhEQWTGQD65lA/ZECeuR91ugY2wR3weiBVYyeV+W19YAS9jBcz2TO+uUdh0RdWEJ1iHTmKNi2PVX1HNW6c5FoSKCpUEqdiZDn81rrVwC01gec9f8ErDKLH5Hku8g50Ei5hBAtgiDytCEtZZFU9dDpLH0GSp3zgB6waqJ40rp45t1bySuWmIWX867jyKECygavTsR8xmOFjL+4krdiU0R9G+gsYhImlAyCnTDz6tXo36ZutUgp5s8ClsF8rfmJ/d9uUKgoDN+Y/N9YzUE0gngI+8nTKQS9aqKk7rmqf6U4P2311h1kIP/KtIRN061iEt9eKJLzJjjyURFcDu8WGzvpoTz2UEJesUjh8aLMiUHaAtVvmpjdg8D9yH89D2iAXUOKmMa/8uDE2+imO6Cd0MInmZZnuVt1bXDpug5xnl2MQLPxwlXAM0CV1voxp/1sZ7Mv4b1SVgJfVkr1UEqdh1hrft9yQw6RDtl453YFBJLnnUAEbpn4c17FC0HpxQkqFoqcGF9bSHxpoTjRxHtQ8cGXiDUU8+4+IZK3PviiR0iHAI7KNA6Y2xeQ8A96A9fB+6u8D5qfOB838028504j/d2hNauYCuOklJclz7KInEt0DkRXAIvwamRGEVXtPLwSZ3k2OfxZLC2Zya6JRTy2YD7xufAHLuAJvkXl+qvY80wJez64MJFBabdfCgVJiv9aHszLS4SNxGPiRNXvnAPihdyOSGRgsvZQS6TAsMo6eh2PM5vnKeUdflL/MazvoF6rWSCqddKUDuURJASmvc+zC9lAG81EpJS6DPj/kcIG/2ua7wS+glhhNPJN/nWt9X6zz12IOvcUovKtyHSM0zkTUXvCPkiX0jXsn0lxkffIz9IpM3mB2VR88CXyBkgoR/xnxpa5GBhFcmq6axAdCojuZJSGvAa4O0/WrUMyD8VACBQSCRQm9xUp7jIY/1IlGxKJXAPwkhL7ZX9gDUJYG0GvAfVtvM9R95ys46ybVH4vkughHx4bMp9aIkSopZYIPTnBGi4X4t+u4FUzdhureIVvTDFk/NeY5Z0B4y6Rdv1c+lOzWKcUnza+WYsikU0K5HpNh/ilcGv+z3n+iDhyndcvxvVqTOf+2Myi/Fm6mrPQipmICpTe1Fhxg6DxrDxNMxFprdcRbNdMWzZXa/0Q8FAzxhWiBbGR3POSdhYkyPMhhJRqhTxBnIC4RRE/ZIjTJcztQBFema6PEJe4sWb+IwUlebJsq6uMMtOqk4ioa0NZEO/bnYq35k/h7EW7qPvgPMhr4NzBMZ5kARH2M3pllUhMB2HXgiKGvV+XcPhRUTMmk0UoEFaFmw8MFXvtQQYmyPMTY8vdwATe2jcB1hry3IQQtf1WPWbOE3Nex/AyDL1qfg1hJuYNOoS6sDFsBWoh7xg8PfSb3Fn6EM9xIwDL9GYv+1NnxD7dKIm2i4duF1PhdqFT7Rxoan3QEA6GIDbB/qCNtPY9fkrdfV7+VnYiJHEMUcNaMt2JRxAxM12HkNl2hHgskVyBSK+AqG3PAaphbRSK+kp4y1ioWz9MjlGcR6/BJ3iVL1FMjB3TRzJn63K2lI5g9PXGAcYmQLgTkUijSCL4iXkUrozDJaBNbltldqkr7UfB8SMAfEIBL3MdEWrZ0TCSkT128NZ6IwGvA1YdlTFyjpfgwX4QWJw6IOvOIZk4rbPUTgjS9qZDDV5K/bZCDVDzpiOJHkSS0gPD5tZxV/RHPJT/A37K9/g7vZnvH1cpDl4hmgGbTL6LICTQDoiUFH1u+r4WTFV2uqBsEgm3taWzZvIwd9Kdz/iM7kKedQiRWQ/aGEKcVtqqQ6RQSyjbEQlzKclEUgfMdY/8EcmohrqouNNt6itEmwc8pal6dTTFM2JsYiyX8zrDS2sYyAF2LSsixnlMOP4Wr+dP4aq9lXC9R5aF6+PyQZAPBwr7AVA7MUIxu6klQgFCoIM4mCDPI68W8dY5RfA0orI9BrAekZY/AgZ5SR66KSFOgFGDYPsBWe9aVMaSvEz291/ZpaRNTdiaqAGi1ciHiFV5fwfyvgsPDP0R8Sjckv80o/M3c31nlUSt120GSbTNNQVdTAINq7F0MmTKhRu0jUVQVpPTwQZqVbdb1o7gVp6kngKqnhktEuM5iBovDyFDK3FaUtwOon4dk5wc3Vf/M+G4AybHrenr2GZzEJdIbWKFPvIzSgkBWRL6DhRdvIsvspYYxVzHy0xmbcJuCSa+cSjUDelH0d4j1A3pxwl6JhIh1HI2BxnECXqxgQnsYCTvvjvO+0g4B5E6XzXnsh2gAiFQO17/uI1keuqAEKmFvVY7zTkMAFZpmKuysoNyg6JiaRbbtSKSbONDkGo3QyEehVj+uRQf38OPe3fij9LG7KEB4S2tZgPtr/SmK5u+n1qS2QaqlPpb4BbE52YbcCNwNvBroBBRBd2gtW7TGKsu9K1weiBT5qFMeB24PMtt3aofHRWJMlcL4LpZS1i+Yo4nKdp6y27MZrFZLkGky96YKirOuhSYCituWpBi9xguCR0ERgcP9hhSueUy+a3rfR6rz+nFdT1eToTJzOaFRFjJ4Yl5HGQQNZyfkKxHsoPuNFDL2WzkUk7Si56c4Hxq2MFILr74bd5dM04kzivwbJ4DwAvDNobXxLwPpwyxbjcSabdBnlr7FN6HxY8V+o7gU03BXiEwt5ZpW8OWhYtOkvFQCcyFF/K/wm6KuTR/IxdZPW9nRBaSaGeGUmow8C1gpNb6pFJqGZIudhqSTvbXSqmngJsRX/U2Q1hQu5OiqSEr8wqz266jk2fZJChbA1v2jWDX2iImzFrD8t/MEUkzjkeeLmybtYHa+XnA3DFCqMecyYanuFmGrMdtDEdw86twXXwq6lHwUv9ZO2teAwU96nm54Tp2MJIL+AO1RBKTxXA+4CCDOMggXuVLvM+FrOZKvsSrCbJ9mlt4991xQp4mm1HCAeiUNvl5X0Okz1jAOG28jyFXV/o85VyvUXiewTtBfZBlnct2JE4Lm5C/+k1DpoXAeihlK/ud693pi1l3lEQKrRPG0g3oqZTqBvQC9iNR0bYa7z/j+Y+3GUIJtBPDb9/IJJG+drhxT9zhzm9HIlGrglu1dgo9j7xE/FgvWJHnEWNQuAV4+WpjCDHuBE4dhbF9hUBfM9sVA1chjjbFfT2esc5GCRwAPoJTo0nEfAKB0hwIAdnxXYFIoADxHjTQHRAP4RjnEWF/Yrd6CkT6NNhNMSfpRYRahhuJszuf8QHDOY8YdbFhoqIehXdNtuOcIEgotiXLgWbMAZn1D9n1yukHxAEpKtftEOjzsyObTxES7ZXV1q0LW1MVACknygQ2cEhVWoV750eQd64pj9YmKurcbaADlFKuHe+XJlsdWut9SqlHkBxbJ4F/R2wv9Vpr+5QFpoxtbYQS6GmEJDK1ql7zu5HU7CXD/VOWUmpbomwZsAB+uPYHXL1iDfGnCmFpnkicO0mjekWIz3rYvma2y0OKXR9CSPM6vNqSA4DefQMKRR9FiLMaT+Lc4t/IaXcznRuptQQh0J2IlPiaou6ZYRx5uoiqfSN5gdk8wu1spTSRWu8zevA+F7KDkZw09LODkexgJH/gAj6jO28wWWI87bAewcvR67fjAkKaMYRM06ibAXE4Ooon0ppzGgWsqoZNoNTmDPt7mKs1xRiVewdAdBIJr9zRleLOPHeUGNVOZ5RH2ijNZ+7lzA5prcc60y8TXSp1FlKk5DzEVz0fKAs4epuL4KEE2gmRzhM36AFJcjqqlQcpXcWHmsPB7e2BsklCms8w0lNtWunKVdO6BLqdZAnsOpK1lnE8UrWetiWI9DYAIdXtR5GPXCDeVwj3lJvuxo+g+mjGk9XiMud4280xB0C/kjqKe8gAbS7eHYyklHeoJUJ3xB/CkqqtoGLVjvVHCuA1JWO30udajHOTO9YTpKqk18kYu11mvHCrZdx164AvAO8lb8cBI41Goe4oTcEYrWGy6lgJPT8G+kvhgFe2lVG+t4LyoTCzvcfVEnBVuUYabTMJFFqDVS4Hdmut/xtAKfUKMAEoUEp1M1Jou6SMDSXQTgLXaahZKfuGBDd3JOnzQl1E5drxbGBCgjz3zC/xAv0HkJwMAbwkATE8O+XP8LSYiVAULV6kIES2Duk34SlqyLPbICHOU9Wm7SQEKvpcojqAyDJOhiI7zsnOZuuA3pojHw3i3XfHsaNhJDGK2cAEAJ7jRg4wkF6cSMpXC5IcYTfFrOFyzusXk2w7x5DXxzX2utgx2Ythl4+SXEX8AJxahyc1W/qw5OmeV8uhPZN6VL8pDk3V18vysMo6SnmHw0PyKP8wdfsOW7i6iWgz8mwdG+iHwDilVC+TWnYqsAN4A/lMBvgqsKIFzyQrhGEspwmSyLWRB/7eQlAXycvEEmc66bMtbaFly6By1nhqibCBCaxmGntWlIi61tjekn7X4jnoHLOkFVAskm2IFGVz10m1FO5Gvm2L45Lz9RyESI8tR6QwEDIxEhiQnEkghlj37DordU40v0biuw6REF/GUxnbNIHbFf2uqKMh3p34sV6MGLyDs6mlcv1VUBzn3MGxpEoxiXWQ7MO0CSHmt0HMQ1ZdC+kJ0J/m4CI86flAmu3cc+2J1mPS9J2MxUoxdw5Umw+V9rSxuwQenYX3n4BUO54O5Y1UWGwOIS02z+fctnj3pqnS0mphLEVKb5rT9P3Uo42GsdwH/CXiWfAOEtIyGC+M5R1gjta6IYdh54xQhXsawyXTlAf+go71VW3JE5DctQuvFZulSxIxRI3oOg1tP4rY7CA9UZxMXiwa5NkKbWzPU3H4WZ4Q87E/x5PCLJn8ufm1r34rxZ3IcFaG1F/u6415LUKm55BIyH5kZ1EiuUNVbDRVjBZSj+WxJ1ZCr4neMSr3XS5jt6rqcxCv20R9UhDyPGHGmMm65yfJA6SGuQx0tnOl7YM0JS3RXK3hBkXUmBDa01HNOhMlcuY+5aych5FG4yweGkByxiEHMicpSEeSR/NU0jYxOnH8aRBaKZGC1vpe4F5f8y7gT1v+aNkjVOGeBshkE02r7n1fE42KJOoiqEwStK7arWwZifwDT3MLFff5yNM6AW0y89uPpqahSyKDk77Jwjj/1Dnb7gRieTLfG5NpqC+eBAqpEhkIc4EQiXsMS6w+SdiS22S8pPAgUmhvLefqJhywYwKq9o30fjflJZOnew2OgUiflnCzcY0ZhJUmPfiJchBpPY2zxGKlWLwU3Pz67anKTZBnwGkVVsY5PCSP6TpPxm0mKhX8DMqXpSfP5UqxXCmGIwr/5UpR4Uy/a4C5TmKH4lY4t3ZHWI2lbRGqcHNHptR+WaX9G6fQ1Z5KN5rJyajZo01F2TtQWTqeepPL9cUVN3lxjOB5k1pbplXVdnMcYk5V4BHcSVJVqi7mkXhrXmXWf99b229sHUfyipDyt8XSOGoMbLceqQcQta0bV+lKY5ZAo4Dx6rWevTYpwWSEREeZ/2W78s5vDl7JNEu6o8z2ruRpaneyHScr0hgk49AJM6ZsbZeWQPuQyigDnbaepk/7wXApWgc5Q2bADSqhxrVoD0k0SQINqhtrYT8o7TPhVMLZfL1xkDJYrhQz7fJklTF5hP+faVF1rlXbZghfaTUVbkTpTV9r+n7q/tO0GkuIjo+0Kf2yeSjf1qhxnlqpTclzjZAnGMlzzbVCBjGEjxLp9pxlP05VkOz04kqD/tfUIERfZ/S2q45CcVTSNNkE6vZYsWKZnzvGlAC7zJDoIEgU6foUIR33DezGVx6FASaudB5eIoU6kmEJcydil7XZkmybe942m5Ilz5jrEetWDWyK488B5KRd8jyKd66WOB0P5WbA/5HW0eKOk1S6LszlqX5Txls2B7hBsXwpzNSamVpTrZTU7FyrIY0vQgp52jqeLeUp254JFbpYLtxQhduJ0VKJopdvBPLlizydCrelUTYLMf0DbzCZip9cKyrMdQhhLD4Kbx8VidNWQjlVAbwCrIdTm2XiPYLVtekwkCTC/QghItfWOtf89h7DvYsV+lOMFGnVnEGSmoX7eqz2hFT7UrZk+BrwiJLJhqFch0eeVrK0MasWMTxpNWaOkaQuzmSTzQQ7bpuFab5vfXBdFaUylvpNQfVSYEhwTHJbogagOiDFoF977U4OhgMVS2XqBQkVbY0z3xhiGMlzn0wtbQtt09AVi9zjQDslOvHQuzZassrCNmDmx8CjwF2ifEwniTYXZZOgcu14ljop1B585mGZiZmGRJyhcQ46Bd4L/gCep6jdIVvYMA7HvmnT1Q0AtkLPiSc4YkJP7j3mvASPgSQfWE8wTuDl29mGl1TewVpEfeuWTgvCWvNrpfHeJKcgTDofPw6maW8MJ5Frmm5fV3WbO6JaSzHx76Sua2tJtMLe41ZpcFAItWyWzCfl0LX75Jia8ABClm1VTald67V2YptmUxESaCdESz4c9oFevhFmplNd5QBXorAvxbIIvLK2jA1MoJStvMBsCc94FS8JQmCQfpCKNpbjyGzYiZWoDsDOQV6ECwhxjRsDb8uL7mx2wRXDEl6zgj6mH6vGtbD92lpaFua8DvX1cvK6qlpwvGgNbOIH1xa61sz742ABT7rOlegOIJLnepI/FNKR50a0LgdEEs3aHjpLw6OqVT/UmgQjxFtCrXgp+f61RJoruR/Au1vbgtTalTy7mAq3C53q6YOWJk8wUugs4DFZzuRMlAl+VVw04tDIz2Dm+tXJ2YPWAavctHAu6Vhp06L5EpD01wfxtO0pv3kkHH3q1gyj6N5d9Lj3M8rRFLOTkw0FYgftDRyz0uvvnPGcNH31co4RlKEILxlEHkKilkBtpqK1eOrk7XhJI2Kb4e0oCWI+FhTvaqXIWFZXInj/V8x8kOdxKpQqBy5t+qFuI1AKbQ9UBNznSXlzaTp5lkVgsek3RtuqUts1LKaLEWhoA+2iCEy2sLp5fTZmx+r5/30sds5X8WIXl4J4j4LQ+HpSJSBoGfIEjxQMwRUjtkebR9aoSvcsLEHti1PfUMCRp4s8O2Rv6/3rtwm647NhH9VIWInzQrMfD0FVYzaRLJGOwksmlBXsa75nxq0yw360pHPMyoysbaKPEpgVqz1DW/ywVVyagrKIR54xjNq2A0Q6tCm6UBhLF/pWCOEnzWV6M1UXSVJxfVhR/gsonxS0Zypce1W6l551FIm/DyX5VcSVmy9wMywGL+i/BiEe++LOjTA9lWK50zoIreej1HJE5WqOVRSFYtDPwcjntlB19WjYBHUPDpPdHszjSF2RSIzFeOEox/oS/GqNkRz2ESPZu7WnSKB1CDna2NZiROp0CDzhqVuEyU8L3rV5j2SSPGnOy71+uSKX674RSOt0Goy3hVSi4xQVG71mW9igI3nlZjOWsknAAlh8fdtLnB0KXUwC7UKnGsLFfaM0I9jC+G2VvHXRFFShZsm+6+D+5Sk5AJr1MntbyHPPihJEGgPPayPm6735jirJxAkeea5D65moM5FQkQGgF5h9ngE2jYaxeMnlXfTGkxyPgWH+NHDDWNx0eqNJSua+HXEosv0GVZVJaxOGZLHUkidkq3ptTVgpNCub6DQoi4o3a2eDJfo+wOI3IWZspV2WPLsgQgI9DZFUgcWWNXPalunNjGALBdSLE891cPO9TzKWTfAL2cbaP4PIM+u6oRPhsSHzqW8oMKVurSR2AiEAG0/ZUupZ60Xq9af1fPMrxTj1/6TupW8GNVl7lU1svt08AA3WG/fYopLC5QAAIABJREFU4izHYc/TV5UFTaLO5lwkHMV647oORImi3tnAXsOWthXngo1Ye2jWDkX3iFduWWWwLbKjomySkKb95AtJ06CLSaChDfQ0hJu+z19NYpnezEh2UEA9I9lB1b6RFN27ixPWAcaxS6UjxxqS7UNB20UnwbPLvsL9Dfdw5LoivCD/mPn9lJYKjbAQsrTOQTIptQgApTK/4Pqdc4B+t9TR75Y6j0TjQO9M5JnJ1uiG3cRIeAYVm2abeOtlvBJkxzY7SfEbg03B5y67x20vWHVuRVb20J/Ye/Pr8hOd1LHUt0Gw5BkD7uiKNs5MaJ1qLB0WIYGexgh6sAuoT9SV3MFIigbXUrdiGC/svJmShXukUFCOcEuixf9Vai3e3uMRx8s2Zn5rzHxLSkpW8pmP1ncAX+Bevo8lOa0zG+hu6fE0xT1i4jC0dp2Xuu/Y4kYkz57OBKmewy62iDetzR5kpc5D4HkfB3nXBsEmtneJtL3JMxnZSKF3aC1hLZ0EZbOAKV3czpkJYSKFEKcjegIr9BpGsoPZvEAtEW7Yt5SbBz/N5BlrwdRHZH124StBjkO2JFq0CnblF3EV/0rVrNEk2wxb/iVvHYeS28oATeqaYOxgJH84coFZMuEidcsRC9enzm8mpJNIbZUTJ9H8Kut5bAtd90W8kTcjeXhBKsC85/QRBL/TUHupb12IKrdJsaGIZJcuUYE1F7SXg1FZBKlSNw/Kp4bkmRFdiFVCCfQ0QVBYim0bqr/CHU/Chp1TeXrnN5ny6FsU8AkAz6y4lb86vFzqIFbhlbIMgPWQhGT17fBCT/qMzoGlJTP5O/6eqvtGw8tWjecmeW8urMR1NTnFIAagFyeIbyo0QvJH0O0y6GYLTF8EXGZ+h+NJfX41aqbzswnerX3USuV9SbZ3uvU1f0dmidbf3hHIM0fck56QXNL0x2e2BcomAUvg8LK8kDwbQxeTQEMCPQ2QSA/m2Dtd8px9/EUOL/DKY/FbuOrRSsYP3sBvZ0xFPQqMMFOVbOL/yvcTp9tec1im4YVAf7ju+HKW3zcHEf9aQ60YVJ8yd7zCNFFrl8QlZ223MXBKwymbZB1EuoziVS2xk39M6bDNN26bcs+SpzsfnHc2+Hybkge47dGkXLkLUskxSOKs8a13p5ZGgjyn5LFSxUPybAyhDTREZ4N1GkpMPmn06fybZWYeMB0R3KKwYf1UrrqrEtYAH4JeA7o2DUlm0abmAbPgofwf+Na01st9IxlF5ixxrckg8ejg70oIy4NAscIjTWOfTCRRcFW1lkTTkZ6LICncnzTC2kDd/lzJ1f5+SkclzVxRfr04EV1RaIiLVLL0o4ZUp7aWQlkEuEfmV6p4zvmduhRCCTREZ0Emta1FLQHlVaqB24HHzfxQmVQEXjucurnf9uQnz6jJnrProSL4GKoY6RR3tsglwXk2eAylMvcdmHXJhwlsoIbhUqPzMiQ7UQJ9obivVGVJqFgPOFMf4KumPVvPXEgmzvWIyjZdIgQ3NWANHc1hKBUbsR652cJKd2qal3/WRa7q2+F4GYLcEKwrCjPsBHAnUAVHo6HTUJMQEmiIzgC/tJmOKGbzPBuYAFPg6IdwdDHoj0XiBOBDEqpbP9JlhrEvoOGFULEdostg2F113D79AZavmGPCNA7iSVItLy1pXW6mvoHrleGnTC++V5jGs8zmZcOY489/w8s8BNxLP7QGvRtUH5DkvTGz9yDoPVMm1gPX4hWnTge7r02XZ2NiTxAcmjIIsb3afdzC4Z0BTQtrYa38V8PTkFsuqfXAc4yLGhKN3gcq6km6SftMgrL7YMuCEZTfCo81hOSZNUIVbojOBn+sp4sz1aMAXLWykqOmaliffPmtOQzVL0G1qQfqh2v3dL/c7cut5rBIrGV3Ilw5Cx794C7Gz6iUlHTdsvfAbAkotRyljqLUUX7InegstLv1FLCBCcQo5oWG2bz17hTJ03sOrNFXUK4103hFKrIkwaQEPLYYjrmp9r5K4yTqSo81eNKlX6p0l1v3Y6Q1kZodqhEs0ahpWUiIWaCiVu7XaAR0HPi6mX8c2AtUeyQLDqHeCCtVVVJMdYgsEKpwQ5xOOImQBEuhbyn07Q8qL1VVW70RqrenVqDwez/a2oiWRK8oBC6FVQumsKp0CuPPf4O31kyRlacg9wLPqbDhKkFhK942faGoLw9e9HDKuh9yJ4/xjaS21Uzj+SOzKWWrNJjaoCNmbGHqmg0AVKy5lrqFw4xa+gukEpkrF/0zqRmIguCS46fOVOPM2747g8q2MXiSaLZQ01rGMUjlARPlg08/DkwxauJamWcKlC0zG08EHoJdQ4pC4swFXYxAO/HQQ2Rj2wP4laqj/ElgCYkQR8jN6UJvM/sar1u1HuJDoJ6zAIQ8D2HIpnnQuhylFiXS8aWuK08iU60l7ETvB/VFaVPzYeaipXxCASe5nAi1bHCS3b7ecDnxnxXy4mU3wTGY+ZdLGaVu4L5DmvH3VlJMBI6VwPfBI7FBeFfPkuh6kkmzJaREf5q+zk6iTcASDXcponfKLVuR+j2UNSpqoWw9RC8FpgH9gTfNchSYCrsmFlG2ro7yy4CHAeooD/kzN3RilWxTERJoF4H+Cqglue9v1WlqHFSvFvqIRoF7YMeyETzMnQCcO3Une35T4hBobqEmnrSZSp5ZYQCoFTB+USXLZ82RtHyrjpJXf4or+61m+fo5iaTyvFwtlVnuhhP0EslDKe4r0eKROxYpuP0yiLq1J2KX/D3wp8BvzUFjCIkG22TTwy+x+gmzc6ls02NjRu1BCh7SsFfBh0KA1W9ml0zBH3JVNsvMPAo8B8Rh8UuJIcG9EDPlb0Kps5noYrlwu9Cpnn4ICllJh/v6Q/mlQHVuSbsTcZ6I122NTUM3Ca5kFXULh5E35zDx0kKn1uXygJ5aDhlfxi8Dc2ADU1AvG5G7KEq8YDnLmYM4AwFFlyXKmgFUrLgWtRS01tw3Hq/k2CEQ8qxBpE0bkvJbkpGN+jaIEGN4ttPTjTg9NDU7EUO05GdeYqTRvckeuunIdHihF6MMQAQp4j0HuARit4dk2SroYgQa2kA7MbIlTxB+qN6YXZo+P+xLSkUQj13MV/2jcMuCn3OyoReshfgVhUI0idJcuSc6SO94kl3mIa2RuFdA6yg/1z9B77dq3nXehnXVUjLs7WqRNl+V1VYFTN0ByVlbV41cBRsyk406NZvYUD86dmKE5qNpoS0JDNHyfy7wYkUhOKlCDR5xDi8EegAfA/dLW/nQ3IYQIgt0MRtoowSqlBqilHpDKVWllPovpdS3TXuhUup3Sqlq83uWaVdKqSeUUjVKqfeUUqNb+yS6GpIyDjXhKzoa8YpcNxXDkSQLSVgiOWSPHCowtStxpE9Ijl1sGtJJl37pJZOHp97vzX9zzdM8yS2GGL9gWi8jUT9zVFRUta8hat21m03R50FIJRk34UFj5OaSa1NItGvYOJWqaNLHXwITTeL5h0A9BWVzTCiVs0kN4kkbnSMkqqYhyUNmwZaSEWyZPiLrw1UolTIxWMF35Denc+gKCMNYknAKuE1rPQIYByxQSo1E3CrWaK2jSC6b75vty/DSt/wNsKjFRx0iJ/VTRW3T1bdu2ArgJf45DlwE01gtGXxs+S/vaCTXxWwasg19yNampqfCrTyNfgOY3BchTyd13nbgHMTWCTB5DBQpvDJsNuPQ78nufNxEC01BVyDRjc1Tnxoi3bWkiN0fF/E5XcSFWn7LlgGzgEr5WCxfCszSxKdCMbup5Wy+rM9lrJ4STIBbFbxkyDIdXgKmQHmOH6OnNUIJNBla6/1a6y1m/lMk5H4wMAPx2cf8XmPmZwC/0oK3gQKl1NktPvIuCCt5urbP1vwKtl/2yn1RjENeUPnAdCjlHW5b+I9iJ8xDSOiUJZ1MsZDtiAGYTEN9kdjN92AUMuUhiRR64xW7JoYowXMhxBCthWHsT0y/UnX8StWJhLpA8//ae//4qMoz//t9bS0EtEAjhRANDWyyBb5YEfMogs9+KdDvIrWoVXys4ha/+nX9tVu32tWqa9mWuu1WXXer1bVYbUXbKrpqqfJsEfPsCxDcaEFTwCaVUSSEFJEgQsJi7+eP+z5zzkxmJjOTmcw5M9f79TqvOXPmzJn7TCbnc67rvn6wyy6eUFcda3iG8xjOYR7jYoZzmKW/gXu5gmXcyFpm8jQLAHjhwvSfGb8BvRWoK+7/nxJ+cpoDFZF64BTsRMZYY6yTzD16/roTsCnKHu+6bcoAyMdlGyTfy348cGiBdZuZhe6FbfBa3WT+mb+1zycBs3FWaC5F1tPjWaFeU+xCYZ4E8wWxxe5rnBt3NLaQ/AVYX0or1o0L2EjbN/CDfIp9Y5DP3Gm0yLm4Qj9kU/BgDfP4DrewgZlsZQqd00Zy3ZMPcdtdttjIlza/ANf3/1ltHdibyJ397VmBVJgFmvXQReQ4bFjl9caYA5L+zivVC31+2SJyJdbFy/jxOqufiVTl+nKJwF16ug0gyiXvc361TVmhDps7dwr8uO7LTGMzHbeP44e3X8tvmGYDiMBabE3AIwew+YsxClH/tlAXWjnJz2EFYInBANIMN+z+BQDLe6+gu77GpuDMB4471a7H1mFvCrwG1sW2QivDyrV/W68RenGrVo3q7QQWMGXoVgCamc3zLODeRdcy4fhOem98Oe6sPwSc3whL26DBnM9iVnKgShjRYy9jj4jQ3mqtCI3kTUGE5zRzJSsBFZGPY8XzMWPM027zHhEZZ4zZ7Vy03oTXu9jLrseJQJ+ZN2PMg8CDAE1NTforTENBXEQ77XxQe5bzn2d5ga512KTznfD0wrMYwhG2MoXreu+le3WNtTZbsIE3y4B7ILFAemHIKXcwDZPfeI0beZ47uS1h+6svTWE69qK6f+i9rNk9j7fPmGQrEh3EunVje7DiOYJKEbfBJuf0ljzwxBNsAY2ZQzdQ+2Fnn/2GAy+0wTerQTY/BW3CiEf7/i9eUtTRRhRNY0lErKn5ELDNGHN34KXn8FtQfAV4NrD9L1007gyg23P1KrnhzXcmbHPPk/+Zbz4ITWYOvzMPJWwfhnU5ZRs8NL/apbtssq5OaoEe+NJ9LxCjnju5ke51NdDUA2f2ONEEzj0AzQCFC7r2isUXwgrd+uSp3LXrxj7bPfEEmMZmYjsn+/XeYwZWHQA+7zaMwc4Ml2uKSSnw01py6h2aA6N6O+neO4pDDOcQw3l510y6t9ewoXcmjx/75ZQdiMCWu3zhFDBX+fOinqt4iTG2ktGFwiMZalFXHBXmws1mDnQWcCkwR0Q2u2UB8F3g8yLShr3CfNft/zzwFvaW/UeQVHxUyZpU4pmyhdlfQdWFcPaTa1nGrYxw1XuGATel6DaRifZ9dr6zodrl0t2BLf23EOqJseX3p9tAm56hcFUVHD3g8j5HYD31r+G7b0MkNB/Ct0+4jWU4EZ3R93u8juU2nhyg9VU4RogHGdGAdd8Wpom3kswYWFpYC9R2gDE2zcqxhzEQq4IYdC+vYRqbU3ZkSThOVertj2wC3oEli1K/XpFUmID2O3RjzDpSz2sCzE2xvwGuHeC4Kp5sg4YuAdr+zeW/VcNEOjl30b8z8QngvsTyZx7p5kLPcn09eQdkLjS+A4wHcwfcV305P+QaWCYuX1JglRtbvMhAzB2gMOIZrHWbXPfW3+d+YE+CpZrS5TsXYtQziv3Ii2A2pvlelxjMEhABjq6zQUad6/ADe8agbtziYL5Z4OOZs5AXgXVVcCJsWTfDuuZbsXP2txmmP7mNF1L0Hg3yQkfq8g9LAt4gnQv1MRU0B6qViEJKsqs23T9o41/56y+sAJ6DCcd30nZhYsmz5IotJD1vwHZjMR3YFBVXcYhqeL+6ih9yDdvWO/fsMmDVHpcnWdzUjmQRTXbnerVyvf3SzZeuqpvDfkYRo57L596LfK+/T3aWZqdXtcibB/0gzf5K/myiGJa9yAFbdGoZNkL8emyJx+24uXth7iJbitFruB3krNrU/UKTUfH0MQIfHZP7ElUiPPTyJ5t/zBf+LfF52wrX/3AxEKgbmsoCDdYRbcfOf7bvsxUw4syFZj7HbJrZtmZ6oEzfWIjHX3hRt8Uj3VxotvOjz3AeM9nAVqbw0H3Xwc2vIj2nZrB6xmCtaK9q0WM5j1nJhU2IZB+NK2IwJrVjTOQR4CSoOdXOy7e6qYXOw4DzKBxzJrTCj/g/THT/K7QFBLMOTVPJB4m2IOZKBZ1qdBioS2j1Pgh07OrTnSJIcjHuhmrsjLdj1aI5DKGXLsa4lA5s7VhG4NeUjQWOUPi5z6Br1rNCky3NbCJ1tzKFx7ovhuvuB+rtxTWtgDbiz38GrU513xaPTdhCZomkEktjBLkMzMPB/V7A/hZP8lOQmr3fqMth6twDNec7z8KZ/I/u33KY4+Mt/qjDlv6rxeYFr7WW6KY8akhXIkbg6MfycWz+seBjGQxUQENIKvFMFTyUqqx6Jletty1odZ41FRo6XN3Qufiu2zrYt7CKKWzl//AjXu1tssV4vC4sHMCKSbLrLUSBQynoaa3GWpdnxTuwpGY9fjPwNzLtqBSQ1HPdq0klrDxi4GFxwvm623gYOAkOer/Pz2K9B4EWcZ1uPns5NvRxIb6A3gpvTarhI7HulcZaOPCefUldtf1jRPjomHxk5UjBxzIYqIBGkKXGZK7VSXrxDD421gL7bJRhY6DTCnXQMxe6GMsa5vFqbxPdy2ucoeldmFKJ58ALJ2Qi6K7NGDCUgot5nA5qeeie67BB4wfg4AgbjesCisbxFp1fn4j5PliR7UItzsEntas+VRhPPYn/Bt7Nm3fDk6613BvAJTAJepqr7Z+5EVhoxfNhLuPb9/4j3OF3L9K6t9nz0ccqJ4pIBTTkJFueXipLdk29/HSUYEF4+SrWkpyKTToaH1i/BdYuPIPHuZjZNPPXzy63QRersY+0YQuqBwlPyorIC3yTBX2shZeYzbJFd8DK+4k3wV45FqkxIAewlYbet1bp98Gez+uE5byUVMTSbPdu8IaR+ibIVR1tApp6uOCcR1k52c5bTNzeybd/9Y/wb8Cx2FxPr1DaScLSVrVEM2EQPqqgUkQqoCElXWJ2NgnbXjAQ2ELwDVgr03jdUs7E1oZqxc53tmETknrhrYU1zNn5MvV1OzibX9n992LD/jtfpS9hEs+nXL/Pvhe4b3MHy+rBXjzXuK17oHMp1kqpB2a5IvPutZCcl5Iv6fqqOlFdPYKahzt4n1FcMOdRVm6+1M57boQ259JtDOZ4vmFYWuwhRxyDcLSCBFTTWEJKf8WxN5HaTQu+xdl4OtZF22hdUbIIet/BXiQAPKu0DtgInA61H3Zyb93lbGAWHb21Nk5oNe5mP1WVodKLjB9YdH76fe4D7nwKGxA0L+nVPdhv9Ndw56uIpLpRUMqDscQt0Pkwli4OM9yK6LRHrSt3hhXOxqCb5wm1OrPlI47JeYkqKqAhJ5suE0G8/p1SC6bNCmfbJmicCtRC1QNJb+jBd1E5q3UU+znCEFvvdi/QesC5b7vcMgZreRZ3zjMb/Pmyevc8Uzm4zwInkT4oqAs4FRthnD5wqBC1eZVSMQY4P95AfbgLFNtNLWt651lvzD3GRuHOcG9R8cwaz4Wb6xJVoiv9Fcow4CYX7NOWoolNfK5zPLS3+lWI2lqhcTw2dLcaXOtDSy9wISyruwGAzUyjhSb/9ZoRNvw/TrjK2QUFLW0e4TKwgu9Z0bHU+9UDsTWpX1MizjA45nzbes+56l9+cQ4jz+yke3sNhxvEiud9Yi3RtkzHUlJRaXOgaoFGhGHA0lvgpsXALVY8k/M6Gwg0v34nTcm+DuKWJr3Yi0Qt8YvFfkZxMY/bO3OvYstesHfuT1PsggnF55/w50CTGQYxWxpQKUc+C+e61TuBc3tgL7Y5wmKo2gc9HwrUwd3XXp2ih5SiJKICGhFuWgRtd9hKQ20r0tezNYF/eq9EX0O1s0xPB24FFtlgITMX+BD2zaqCRlsr9hDDuZMb2fbsdL9OwtFUn5QuQGPwyaoaUedTWOtzLOnHHUPFs5zZY+fzJ7mnz7gq8ddby4k6Q9WLNoULsDeW6r7NmWK4cEVklIisFJHtIrJNRM4QkWoR+bWItLnHTw7C6SWgAlpC+qSoZIiwbXsy87E8ofSib4NCGi+SAJhjYe2iM5i4sxPZBmtvP4Pn+QJX1P0AsCL68otzbJK5Z33WgO2y4lH6uc9ksivp9xMKIZCFavKtDDbub78aONhmHy8yVjzfMDYneCdUuY48T99T3P6k5YgXhZvrkgX/Aqw2xkwCTga2ATcDLxpjGrF9lG4u2omlQQW0hATbkyU/BmnLInWlodqW8PNE1IvEbah2ofinAwtBPoQOrJ/36VlnsZlpNDObDmrZyhRe2LXAWp4xrHh2uiVutYVPPD1SFZtP5CTCYjUrJeJgIMhslYGpgf+tD4FH7er7rvEA12ufz1ywc6CFjcIVkRHAn2P7UmOMOWKM2Q+cg70rxj2em/oIxUMFNCT06f0ZEM3GF/3tjRm6Q3hzoO0kFk7gdOwc54ewtu4M3uTPuKbuLtYwj3aXDFNLBwt4nqrjDvmW52jwo22j04XEE9HEiNx2bAGIVJVplMphlhPRdlgift7vSc4KBf5+4TdKNrpyoAgu3InAH4CHReQ3IrJcRI4FxhpjdgO4xzHFPbO+qICGiKUp+gt6Qtr4PxPFM5gD2oBvcb7gatV6Vqi4eU7mwqpZc+iglvcYbe+ugYe7l7CBmXRQS4x6elZW20IL84GrgLNHYNM6oAS/z7yxeaFBF9xYit16TQk7h7GBcK8TD4abhxVOJ57c4u99wy9+aFfUCs2aAaSxjBaRlsByZeCwx2ADGO43xpyCvaINurs2FWJCUJaqqanJtLS0lHoYJSfTHOjNB6HqC/7zYKNsz23bONXOfbbvC1RQqQYuhNfmTI6LZox62mkgRj0bemcC0N1SY+c9R2OtzxagtQ3fYmvDz42M4QtS+An2E1UUG9M+D5aeilnpW55BtlPPpMhHnKdHRF41xjT1v2duTG4aZn7akq7ES3pOk9a04xGRGmCjMabePf+/sQLaAMw2xuwWkXFAszHmM3kPPg80DzSCJDfK9kTUE8+GavwqQ7UkNPj0RLSZ2ezormdo1RErnl7E7QpgMbZ4QpxU84bREE9Q4VSS+Tvrvm3GNo9PwSRiCY0GlGyRglcWMsZ0ishOEfmMMeZNbOHRrW75CranzleAZwv6wVmgLtwQkani0BeO/WW/7/dyQCXYOeIy6KwbSYx6YtQzhCNsYCY7uusBZ3nuxYb2e06AFWAF0qtz+y52/lBdoErUqQfWW0/L6Ay7qXjmRRErEf018JiIvA5MA+7ACufnRaQN22Lpu0U7sTSoBRoRPkdzH8sTEgvH49y4zMJaneOBd6CjrpYNzKSWDr7e/U92nhP8fDiv3i24ikNj3SfE8DuveBapRrEqEeeZM+E2bHWqf0yzj4pn3hSjEpExZjOQysU7N8W2QUMt0BCRaQ70qNxF41f7bg+mq3huXNa7F4cC78Ep+7Yxjc3c2r3MNpTe7l6v7/Gtz3qg1UDNWCBTPVlFiSr19uEe7G/9Knht42S7bYbYRRkQWgtXKQnZtClb+i/06QMaLxzvCidILdYCBTjePuyoruGCD5/i0taVVjDPBJ4B1lVZ8fRK9tEFncG5z3psxKKHWp9KhDluif1JXw8cJ3AFTJ+xzXYi2mhUQAuAtjNTSkK+TXo98WzfB/IiVjyvci++Zx8mbu60K+sCy178bSu85+9iXbd7sCkrnjs3U/k7RYkQo7Hu2/nYHrdKwamkdmbRHXkFsvSr8MK/2PUG3HzneBCg8VvAO7DxgZMZwx4mfuhE8z3gWPjSsU/5kbbJHD2AFUuvy8pJRTsHRSkW6drMxaOwDz4CzQ0w9Ux707gsaUed9xwwldaNRQU0JGRy4f6NqaJ6Rg9tAfEE16IMrNU5De6edTV/u+9+ZB02bXMudM4ayeNczEvds+2bWrBT8aPdeiv4EbfeulcwweudqdanEm7Si+f9cMxSOHo/HLPEbjwTeMBg9v2JimaBUQFVQsM3XYTte9LjeWMTKhA1LsLme14KbIRrpjnxBBgPZjx8i9tpZjY9zdV+qkoLdr31ALbK0Fj61rj1iiYcRlNXlDCTqcG5MVcj9wFcDdcBs4EHgAsEFql4FoNKmgNVAS0BwVJ96biEQHoK0BjI7TQ9bsUrlvAkMAeqHsFmSB0Pb02q4RnO4/5d18CKKuu+9fI8O4HOdfj1bT2Lcw+++/ZwYJuiRJjr2ognfjUDzIKV6wHttlJoTBEKKYQZDSIqAcEuLAnbvPq3xtj5zQBtHUCd/1xmYBtjz8AKaRswC3qmwWuTJvMwl/FT/hJaUognQfEEa22+EVj3UPFUwk//labagU1ugXiel0bdKgNEBXSQWSqSsQ9o3DKd5r/eWAuNp0PbJuCLWLGsg7eeqIF3sOXIGuG1aZPjvQy3MYU9jPEjcsGlqLThi2eqykKH8asOKUo0yL5c4yYSxVQpJJoHqhSVhI4raQKHloqwdHHSxp1WRDnJ9vTkKpi4ttNWGzoWOB6mr9/GWwtr+CHX8D6j6Pz6RBswFK/T3xY4YFAg1V2rlCdaB3nwibIg5opaoCUiOAeaaj7UPO+vt7kiCaYNeB74DvAc1vqchu3F/h5wPEzY18ko9rP2vrPtm7xyfXvBpql04QvlHmy5Pq1xq0QfFcvS4xVSyHWJKmqBhoDkoKKlxtAmQuNiK6Tt+4A6kPHY+c4qbPrJZLh72tX87R33s6n6ZGY8uYXvLLqBx7jY1rb1GmODy/WERPHU9BSlvFARLS2VFkRUOWcaMlLNg7aaR3lq/WKMK8XXtsI+Nk7Fpqo8GnjD/wTega89dz93L7yaUeznXxf9DT97w1TtAAAgAElEQVR79n/bGKEzgVUGjhEnnl6BXBVORVGKh7pwlaKSbu5zAc/z3qxhsN5an423uBdmYYNjx7vn1Vgx3QmcYjftZ5QVz9uwVVbucfsePQDHjXBPNqHiqShKsdAgImVQCKayeOv/+76fWWGc61ufAPwSG32b3CT7Qvv0Yh5n3H37rVu31bM6DbbIXxsc7EJRFLBFQ9w0hlYhKjhaiUgpOsF8z/i2v8KK5zTgOdtlBUCq3A5V+AIKcKHfKHt+72o73/mIAbrg6Ltup6Bwati+oljxrC/1IMqaKAcF5Yq6cItENu3JAFjv+hD+EuuS3QnMcm3JPLz2ZMcCl7kFGEIvn+teS/c9NX47Mtrwo23BRtnG8j+RCiVTeTglyowt9QDKGi+IqFK6sfQroCLyYxHpEpHWwLalIrJLRDa7ZUHgtW+ISLuIvCkif1GsgYeBTCKZstpQ4PnTLIALBbMQ27ezGuumrQP2+S3KAJuusgiYAU/XncWNdd/mgrpHOf6hw/RMqrbznq0GW8/2Daxgvu6WA/Stc6v0h0ZzlhtfJEE8j1tSqoGUNToH2pdHgHuBnyZt/2djzJ3BDSIyBbgI+B/Ymbo1IvJnxpiPCjDWUJFNPdt0InoJ8CVXCF5m2O2mA2Q98UAhmQENG7GC6rZtnHUyMeo5xHDW9M6zTbH34uY724E1aUaigUNKpfMKcJp7VIpJlAUxV/oVUGPMf4pIfZbHOwf4uTGmF9ghIu3YX+3LeY8wRGTtlu2Hxmexrtr3sLcZTwZctl6921obAuS1YemZBndyI4cYTox6ulfU2CIJMaBVSBTPoMWp4qkoFieexy3BfJB5TyU/vEIKUUREZgB3AEOB7xtjnunvPQOZA71ORF53Lt5Pum0nYKXB4123LZJkKvie9TFSBAzF8UpyzsJamYuwggrgisebNuB0u2kMXfyGaWxbP91u2A4cBCuY9Sk+XcVTKW+MWRpfMrMHey8/zDbWVopClOZARaQmadPXgIXAfODb2RwjXwG9H/hTbMzobuAub0wp9k2pOCJypYi0iEjLH/7whzyHUVzSiWWwnm22x/H2fdf8wG78/7DWZh22o8oirGiOx4qouw2RRnjrhhq+duxdbGAmnesnWuHc7g4e8z7lkHv0rE8VT6W8SRbNoJimF9ST0mxXCkWE5kAfEJG/F4nnOuwHLgb+H7IMHMlLQI0xe4wxHxlj/gj8CHtrB9biDDTd4kSsLKQ6xoPGmCZjTNOnPvWpfIZRUnKxRL350gZzPuPosOLYiN+OrBE6l4yEucDGwBt3Agugi7E0M5st692E6XZsqb5VBngBeAw7B6riqSgeiUL6xYTXRL6HyPcGfUzlTpSCiIwx5wKbgVUicilwPfBHYDhwbjbHyEtARWRc4Ol52BR+sCXOLxKRoSIyASsTZTtrnypIKO2+IjTRwmf4HW9d6zwH1cCHwHioWd9tLU/PhbsIeAC4DH7KX7Kju94K5zPYx1aDbUu2B53zVCqN3NOMgpehLwFXo/8rxSEqAgpgjPkl8BfAKOBp4E1jzL8aY7Jyi2aTxvIzbBDQZ0TkXRG5HPgnEXlDRF4HPgf8rRvMb4EngK1YG+naqEbgZkpByfW9/7DccIy5gUn3vc3izU8x8cJOW9t2rn19X10Vq2bNse7by+C12yezccnJdNaNZEHdU2xmGj3Lqq14grM8vZzP4EVALwhK+ZOLePZJR1pyKpw9Au9/RdOVKhcRWSgi64C1WCPwIuA8EfmZiPxpNsfIJgr3yyk2P5Rh/+9gG25FmmQXbULVoEAKS8pUlcXAI8Kw896jZ2U1v7x8Lo9zCSuuPZ/FVz0Fi2HVpDnUTtrN9J3bqN7Zw8y6DdxbdzkX8xgA7TSwgZm81D2bnpXVNl0lhr0tOUYC1YY8VDyV8mXghS3GwOxT4QqYM2sVa3fNg8VLobkNkWyCkJRsiFgU7jLgDGAY8Lwx5jTgayLSiNWwi/o7QHRLQAwiqdqNBV/rQxuMWr6boRzh1stv5O/4PjdyJ/NYA7dCTzXM+3AtHccmBoHtYQwbmAkQf+x5oNoWh9+Oy/cEv8qQ96jiqZQv+Yqbb10Og3vPhxPh5FkbWbtrHnNOWMNazsYWG0FFtEBErJ1ZN1YkhxGoe2qMaSML8QQV0H5JVzDB254qIrdtE6wdOod6dgC2y0qMer7F7dxY9326GMsY9jBxZycAd9ddTQe11BOjlt2AFdAt/zDDzoMCHPX6eL6Gnfcci98MW1HKl0xu1v5Fb5h9uO57wJfYwi7geNbSCzyFrdylpRsLSYQKKZwHfBn4b2z0bc6ogOZJJiu08Qlg5zb21VURYwIdjKODWhbwK57hPGbTzGZOYX9djGZmAxCjntm8RDOzWc4VbPvidOuuHQ10GmyAcxe+YMaKfYqKElJyqWcb9M48lvSaimehiVI3FmPMXuAHAzmGFpPPgqwjbU+HpU+kf/13fIYOalnOFYAVzQ3MpINaJrOVZj7Ht3pv9wslTMKJZxc2I8gTT3XZKgqchsj9ebzvDTzxVAqLNwea6xJV1AIdIPGenqcDN2DzOt8D3oEjdUMBOMJQtjKFBtqZzUtxN+3jzmsQc1WE1vTOo7ulxs53rnIW7jECjIWjj7hPPIy6bZXKxpvCeAW4GpGnMOb8vI6k1mfhidAc6ICpnDPNEy/KNlPx+Lh4NtoenTXvdcOxMHZfN1TDEYbQQDt/xpv8js/QxKvU7Oxmed1wZrKBDmp5ngX09gzx8zy9Oe2j75IYLKTiqSi+iN4PF6QW0XTi6M2pqngWnii5cAuBunCzIFWgUNzyNMYWMpwBndNG0uEqIWyf9mnkOajZ2c2f8Sb1xDjCUNppYA3zuKbuLroYExfPbeun23QVr0RfzVjsvKcnmuq2VZRE9thlZRvw2azflV3tXCUfolSJqBBUnIAWqqNK/FgzhBWzzueaOlsOePrabTw97SyG0MvdS67GHAv7GcV9XEOMemLU08UYuhjDU7su4DEuZtv3plvLcx02NmiVgc497kk7vniq9akofbHBQSJPlXgcCkSrEtFAqTwBzbGbSr/HuQq6GEMD7XRQy745VaxhHhuYxcU8zo7qGh7nEj7D7+igluEcigcR0VJF5/cm2jzP0dhiCQA1gi3PN7wgY1WU8ucxNDCo9FRaEFHFCSjkZ4WmdOM+DGahbTPWTgOPczGPcBkNtNPFGJZzBV2MZTiHmMmG+DKEI7YhtueuXYa/vh0npMmF4dX6VJT+0NJ8pSVK7cwKQWUKaI7tyIJ8sxrM8YLZJ9AI8mFiNK0nk7NpBmwTbLBWaox63mcUDbTTfVuNtTonuQNPwo+8PQq2Dr+KpqLkisjShEUZXCrJhRtd6S8A+TTHNscLtMGO6homTO7kN9WTiVHPS92zmTAyxhCOMI813MEtbGUKU9jKt3pvp3ZoBx29tTZN5bvY6c0V2EpDXn/PGifonQeA9YFPViFVlHzRMn2DR6VF4Za1gGZKPcn3OP+wDzgelq7rhPFwR/Ut7Oiup6e1mm1UM7Kpk7taboPtMPny19jKFLpX1NCNszivcgdtco/LsRWHAI624btuX0cjb5VkkoUgF3EIuzVmTPEsRk1dGRwiVkx+wJS1gGZDQjpKun1SFZJ/RGA8HGI4Q6uOMLSpk+6WGmthAkzCryi03T5nHa6yELDKHXC7wNFX3RMv31PTVpTwk29bMU8okwUz081A2MVf8YnynGaulO2ZJrQcS2GJBovB91coIfl9AEvXwdq6M5jHGmJD6+noreXkWRvj+21ZP8OuPIMVz73APd6rr2HFcoyb74RE8VSXbaUTFJKg2CSLzECtz1TvL4ZYJX+O9zzb8RfTOlUKR6W5cMsmiChVkYP481TdVHIMJPJEFuBvTBUAwznEbJqpJ8a8oWsAaKCdWG+9fZNnee7FT1HhtcBRPdGMYUXzddT6VNKJirc9WUiKLSzBwgPeeilcoYX4XBVhpZCUj4AmWZu5WJPZHt9bqnf2sG9WFfXYbirz8MWznYbEN67Diud2Ai3JugI7vILfZeUwKp6VTdDSytZqzFdU+rM+C/U5xSDTTUaYxllpVFolorJz4RaqUEJ/bGBmvCj8cq5gCEdooJ0t62cwsqnTum73BpZVyeOKJT1X4VQswTnCYhL8DG89zKLpkWlMqW46wngO5YwGEZUJhYrCTea1usnx9UMMZwpbWdM7jxj1jGzqtDmennC2EGhJ5rlvY4GjHUDFU8kkXMnr2YprLmKiIqMUAq+QQqVQ1meaMbLWuW+HkZt87WMYHcwkxgRiTGA/o+Ju22AELuuw4jka6AweIRZYV/FULOnSLFJtzyXwppzJFFhU7uceVooZRCQiH8NeVXcZY84WkQnAz4FqrHVyqTHmSFE+PA3lMwfaz3xmPNLWLSdhxRP3mMucaYwJ1NIRLw6fgFeSbzUw39voWZ+xwI4qnkpfNMhl4Kh4lpYizoF+FdgWeP494J+NMY3A+8DlBT6VfikfAQ2moyRF5CaL60nYHif5MpMNbGBmXDzjUbce2/HFs9WQGHkLKp6KRzCqNdWFX8Ugd/QmpHQUK4hIRE4EvoAtPYOICDAHWOl2+QlwbpFOKy1l5cLNlK4Sfy7C+Z7YPgtvLazhbH7Ff7CKFzk7/bFF+Kw5i9uYRxdjEqJt64fGbN6nV5LPs0JjkBhxC5rjqXgkB/AEt6XaT1HCjqFoQUT3AH8HfMI9Px7Yb4zxMunfBU4oxgdnomws0GyJp7sYA+Ohi7Fs2zWFtesTxfPvuSXxfff662t657Fl/Qy/WAIk5npOwi8Sz7vuMYaKp+KRz5ym0hf97sJG3t1YRotIS2C5Mn5EkbOBLmPMqwkf1JfBScEIUHECGuSKaT+gnQZuO+Fb8AzI1/3Xvs0diTtfa+LNsLtbauJW5pb1M9jy7IxAoQRHvC2ZoiSSXGVIUcqFAbhw9xpjmgLLg4HDzgIWikgMGzQ0B2uRjhIRz4t6ItAxiKcKVLiALuc6NjCTH/T+jbUYt4N80S7JXMG93LD+h/yqe4HN8/QszOQ+nkHR7Bz0GyIlAvRXKEHJHv0Ow0eh50CNMd8wxpxojKkHLgLWGmMuAV4CLnC7fQV4tpjnlYqKFlCAH/I1Vg+dz22XO5etE8agiC7gaR5afx0jmzoZWpUhSnpv4HE12PnPLuAQfpcVRdHuIMVEv9PS4XVjyXXJk5uAr4lIO3ZO9KGCnUiWlFUQUb7MYDM/5W6u/uXd3P/Q1+Lb5SE4+fKNHHHRtvE8T4+9wPweaKlK3N4SfPIJrIAqil7clfKm2IUUjDHNQLNbfws4rWgflgUqoI4fYoXz/vk9NgXlqiqYBEcYwo7uet9VOxpbJAGc29aJZ/J851FDotWpaStKZrTx88DR76/0RLm2ba6ogCZhTrCCKOfahtgX8zgvjZzNWs72i8JvD7xhO1ZUVyQfqR0/6nYgWadKpaAX/9xINf+pNyGlRduZVTAS9KBPgm0PTWcz09jQPdMXz71YwQSbtuuCj1jsth01zvoMMraYw1ZCTrBIgga9KOWMQfjojx/LeYkqaoEmkSCi223OZ8+yat9F681vzscXVfAt0GPE1b/1rM8PUPdtZZKqQHwwAlctJaXsMHD0aHQFMVcqWkDlITCX97U846I4CdtZxcMrFJUqv9OzQA+61zs9q/MDoJ5SFFFQq6c0eEKZSiQzVRxSskd/0+HEGOGjo5UjKxVxpvJ1MN93hRK2k1AxUR6Cqgv2AdDTWm03eiLqVRRah++2bcEKZFBove2tYAOHgnOfMQZTPLNpyKwUh3TVhdTaLBz9/ab1e1YGk1AI6G8DaR5TeI2tTGdY9z4Oj6zO+5jyRTC/9EVTvg6cC1VT9zFhZIwpbI3v+6vuBQBMnvVavLdn9/Ya29uzpcYvmrAXv0h8UEQPBj/ZE8vX6VsHt/ioaJYOze0sLvrbDj/WAlUX7qDSc3A4ssvLp5yObAdGVyPBaj9u3czq/3heSb6EikKT4OpZdwPQQDvtNNDFGOqJMWFkjB3d9ezormdb6/R4pG33bTW+5Rl023rr3rgSInDHYsXzA/dc699WOqmKxiuZyVcs9TsuMQYV0CAi8mPAK+Y71W2rBn6BndyLARcaY953LWb+BViArR6wxBiT3MurL53YvMsr8HMskwqznzxrI020cAGj4m/z2oqtZLFfx9Zzq3qW4mqSirv7jKGLrUxh264pTD5hK9uene5H2m7HumUXB97gjckTzMUkiafnvm3HBg6peFYamS7eemHvn0JYmcFj6Hc+uBgjHP1vFdAgjwD3Aj8NbLsZeNEY810Rudk9vwk4C2h0y+nA/e4xMx9iZfgCrCiNJmGecmRTJ/NYw35GMYr9AFzBctpp4Fye4Qr2c9b3awF4qXu2ncvc7N48H6oW2znOZmYzm+Z4K7JmZrNt1xT/g9aRnqAF6olqvVtWBnfcg4pnZaIX64FRDBetWqSDjfDHj0Lh2BwU+j1TY8x/ikh90uZzgNlu/SfY0ko3ue0/NcYYYKOIjBKRccaY3Rk/5Ciu8XSXlWuA0WNhcQ8Ao4bu5xDDmU0zXYxhJhsYExCoZmZTSwfj6KB35BDW1s8DEsvrTRgZA+Dh7iUA9DQnzq9uu2q6XfHmNVuAJvdiUDxXY63bg9igoZXAVLfOHuCNjKeqlC96se6fUsxj6tz0IGIAdeH2y1hPFI0xu0VkjNt+ArAzsJ/X5DSzgPLfpAy4aani0+ds51Jn/O5nFLNtGUT280n2cjwA41wXm987y3LyCVvhBLt9N7Vs2zKdjklH6F7tUlKmGt9C9TjRPR4MrI+mr9sWrDZXAQ+4560AT6U+B6ViqNQLdDY3DqW+uajUv82gY0QFdABk3eTUNUx1TVNr+u6wDtgLbzOJDefM5BIeZwx76GAcv+MzbGAms2mmnh285yJ9RrGf/Yyinhgd1PIeo+28JtB93FgrnK1ilxi+m/hdEsXS+/xY4PloXHuyD+CBESSmq4zF+nI9AVX3baVRyYFC2eYbaxRtBWCAo6lkoDzJV0D3eK5ZERmHrxzvAnWB/dI2OXUNUx8EEJnSV2Rd9GvNOW9xmOE8xsVcwuPscJ1RRrGfzUxjM9PY7wKLvLlNjy1bZtiVqe7wrQJNPe7Vqr5Nr48jUTTB7/F59ADUjABG2O2dQZFsA15JdZpKBZAu/7PSSJX7qqJZgRwt9QAGj3wF9DlsA9PvktjI9DngOhH5OTZ4qLvf+c8477pH5z91wtb50EQ659fy6RNicfEE6559kXkAvLxrJief4Ptkh3OIl7fMYeSkTmt5ghVPgBWBudFg0NAkfPFMrjQ0Ghsp3BnsrjIWa2kmz3uq9akoyRWXVEgrBENFCajYeJ8MO4j8DBswNBqrDt8EngGeAMYD7wCLjDH7XBrLvdgwm0PAZcaYllTHTfyMKQbucs/cdOrUU61XtAXrOp0q0AS3PWwbXy/nCoZyhF6GxI/T+exEPn1OsDwQvL1+UmLN2qDVGQtsO2psHVuPo3uAA1DT6IQzKIxekYQDaJ1bpZKtzoFSDGHVv0f/iMirxpim/vfM8bhTmgwr+r3k9+XU4oyn2GQThfvlNC/NTbGvAa7NfRheHsswrIA68RyNleIWsVG6e4Vlv/8On/7TN5nFBtppYDiH2EE9nS9OhOPg7V31iYfuwUbOgp1qjbl1zzULfveUo04kjxlLvJdnZxuJhREOo8JZ2aSb76zE+c+Bkvz95fqeZNTSLTEGGxNaIYQwYcdZdi0j/GIITUCT63LygPB20yT2nzuK3p4h9MSqOf/kFbTMPcLb903qWzAhaHFuT3ot7pL1rEsXHHT03aQdVTyV1ARFU8Wz9OjfoMQY4KNSD2LwCKGAAqzH1mRwHCShQwoHobuqhqrZ+yAG7Sfb4KFPX+sr5Nu76m1pwL34eZpBAe3jlvXmM9uBk5LG4wmmiqfiE+y6EtymFBe19ENOBc2BhlRAHV4ZPrDW53b8qj/roOeCaphma9u2uKoHb7/o3rAdP58zhmsxhv9mIFFAPwise0FBQbFU8VR8vIt4MNJUL+rFQaN5I0SFBRH9SakHkEhAoDyx87qeeMzGWpQAS4HNfvrK289OiheCj78XrAXbDFYE+xNPAq8fCCwqnpVOqnSVoHjqRX5gpLoBSeUe1xuVEOMJaK5LRAmZBVrvr9ZgrU6vG0rQ/drqHqf6mxKCh1pJdN32qXGbSTxTpaGoeFY66eY59cI++KgLN8RUmAUaMgH1mO4Lp8e5+NZlDb5V2gqx3np4psoXTrBppSvb3JPkuU7ITjhBxVNRSkc6qz6d21zno5XBJEQCWt93U7KIgp0DXYHfi3M+tum1J57bgVYvQOj1NJ8VFE8tfKBkRi/Eg0euaS2ZXOhqqZaACrNAQzYH6jHMXw2K6GjgNrc+CZgKZ5yz1j5fR0A8Xye1eH5AbuKp1qfSl0xWkc6Dlob+LFVlEKmgOdCQCOgQbAGFMUAD8Xqzo/GbpgXzOadie3KeCC///nN+p5TW5NSUdHgpK5lQ8VR8MrkGg1aQWjyFI9N8sxJSvEIKuS4RJSQCegRbQGEEthc3fvBPa4rdPTfuauBOseI6GlfsHRLnOb2lnUThPNzPoiiWXOfhlMKRi0s3UxSvMkh4hRRyXSJKiOZAY8AnsAJ3pr85GH27HZveMpXEdJWE6kON2LSTLnx3bbK1qQKp5EamYBWl+CTPcaYSVr2hCQEVNgcaIgEFK3ifwJqen7ViOXpEYneUGqxVOhWXK9oGnWMDO6wHZmEFtIHETimg4qlkS7oLsV6oB5f+SiXq3yFEqICWkj1YAfXWG6G1DRhr3bOeYNLoXLsvuH09qzN4nAbU8lQKiVqditIPKqClIDiL3I4vohDvhBIvwzcWfz7Te+4J6XRsIJInpu1oFSElH/qbe1OrR1HSoAJaCvZgxRD3+FmsOzZVwQNPEOsD28AX0pOwrlvN8VRyR8VTUfJELdBS4gmeZ0HG8HNCk1uJjcVv7glwGvBK4FHFUxkY6rJVlBxRAQ0Dh4FNebzvNPeo4qkMjHTiqdanomRAG2pHhaDL1xPOV0o0FqVS0AhcRclAhTXUDkkhhYHyCr54qvWpFBdtXaYoGaigUn4RtkAhWSw12V0pBP39brRIuaKkocLmQMvCAtWLmTIYaK1bRemHCmuoHXkB1QuaMljo/Kei9IMWk48OwQtZqnkpvdApudDf70XnPhVFCRJJAdWamEoxKOe5z+SC64pSFCqsG0ukBLQ/i7O/9yjKQIjabym5wXe6dUUpKBU0BxqZKNz+OjIoipIaFUtl0NAo3HCSfDcN+YupirCSK2H/zST/b2QSTXXnKkVDg4jKm0yWbC4XybBfUJXCEYW/dX8F8Pt7j6IUBJ0DDTepLNGB4OX25eMi1rzA8iZqf9t8/y/UElUKhuaBRgNPSPtzR+WTu5csqsFtihJWckmzSf5tK0rBKLCAikidiLwkIttE5Lci8lW3vVpEfi0ibe7xk8U6pXREVkCTyXQRSGddZnvhyFZIUwmvEj2iKi65WJJqdSpFoThzoEeBG4wxk4EZwLUiMgW4GXjRGNMIvOieDyqRicJNRzbJ78U8fvJ+emFSSoH+7pRQUIRuLMaY3cBut/6BiGwDTgDOAWa73X4CNAM3FfbTMxMJAc1XBAcr8V3Fs7yIYsEErZKkhIL801hGi0hL4PmDxpgHk3cSkXrgFGzD6LFOXDHG7BaRMXl98gCIhID2R/DCEVX3mxJOoiKmKp5KKMhfQPcaY5oy7SAixwFPAdcbYw6ISF4fVEgiIaD95X+G5QKnlmh5oeKpKDnizYEWGBH5OFY8HzPGPO027xGRcc76HAd0Ff6TMxPJIKJCp7IUkrCOS8kezx0aBfFUlNBR4DxQsabmQ8A2Y8zdgZeeA77i1r8CPFuoU8iWSAqoR9jEKmzjUfJDxVNR8qQ4eaCzgEuBOSKy2S0LgO8CnxeRNuDz7vmgMiAXrojEgA+w9xBHjTFNIlIN/AKoB2LAhcaY9wc2zGigLtzyodxFtJzPTSkhRaiFa4xZB6Sb8Jxb2E/LjUJYoJ8zxkwLTAAPWm6OXgSUYhC1XN4ojVUpc7QW7oA5B5uTg3s8twifEWqidgFWFEVRcmegAmqA/xCRV0XkSrctITcHGPTcnLCgIqoMFrnctOnvUikaWkw+J2YZY6YDZ2HLK/15tm8UkStFpMUmzx7K68PDPtcY9vEp6dG/naLkiRaTzw5jTId77AL+HTgNl5sDkCk3xxjzoDGmyc6dDs/y85ZmfK4ohSCq4hnVcStlRIV1Y8k7CldEjgX+xNUmPBb4X8C38HNzvkuRcnNUOBUld/T/Rik6RSqkEFYGksYyFvh3V07pGOBxY8xqEfkv4AkRuRx4B1g08GFmboStKIqihIAiFJMPM3kLqDHmLeDkFNvfowi5Oenq3SpKofEqEZVTLmi5nIcScoqQBxpmIlmJKMyl/DzCPj4lPfk0YQ87+ntUBg2dA40GYb3Q6cUq+oTtNzVQyu18lJBSYXOgkbRAkwmTYIVpLEruRFloojx2pUyosDzQSFugQUo9X6XCGX2iHKimvz8lFOgcaHQp1UVEL15KKenv9xfFGwIlomgeqJItKpzlRVjn1DOR7jeonYGUklBhc6BlJaBe+oG3Xgj0AlRZlIN4BonS+ShlQoTnNHOlLFy4XiHt4AWlEMKn4qmEmUziqL9dpWSYPJaIUhYWaLqLRTaBRanemyzGSmUQJWstm99nqQPrFKXcKQsLNBOZii7kul0pX6ImNP2NN2rnoyhRpCws0GxQUVTKCZ3/VJTSU/YWqKJUIlEod6koUUcFVKlIomyd5SKMKqLK4OLlseS6RJOKceEqikdykFiUxVRRwkVllSJSAVUURVEKRGVVUlAXrlIxJFfn8fKHo0Q+Lll14yqDR2XV8lMBVSqCVIXiKymNqRzPSQkjlTUHqpZzHj4AAAgPSURBVAKqVATpCmakei2sVmkliGAlnGN5U1kCqnOgSkUQFMVkgUxllYZVRAdCOVQmKue/T/kQXZdsrqiAKoqjHASmnEj+e6Sqda1/r7BRWUFEYkzpK/mK1Bq4stTDUMqU5ItwFC+6hXRthvX8cznHsJ5DVBCRV40xTYU/7lQDT+bxzilFGU+xUQtUKWuC85xRjLotBlG9iYjimCuPyrJANYhIqQj04htONGio3NA0FkUpG8qlJmw53gB4f5dUf6NyPN/KoLKicFVAlbKnHC7Ghb4JKPVNRX8u9VKPT8mXyrJAdQ5UKVvKQTg9yqXJe6pI2lQEz7ec/o7lT2XNgaqAKmVHOV5wiyGegx1MlK6YRbpzK8e/Y/lTWcXk1YWrlA3lHGUb9fPqr2xi8vmVg7VdmegcqKJEiqiLS7aU83mmEkwVUSXsqIAqkcdzRQajOpXwoH+PSkODiBQlMnjimarjilI6VDgrEQ0iUpTIUEmiWWhBKtZ3NpBxVsLfsbypLAFVF64SedTSyY+wfG/lHPxVeWgeqKJEjqjWd82FqOeCJrvaU72uRJ3KskBVQCNMqgtOlC+wuVJJ7ttiUMjvLTmQK9VnVcJNjqJ5oEpEqUTxrDQKLXqFJNPxVDwrhcrKA1ULNML0d9dfrlTyhbgQf+tCf3+ZKgxV8t+qMlELtCCIyHwReVNE2kXk5mJ9jlJZVPoFeaABN4Us3p6p002l3dQpHsWxQMOqJ0WxQEXkY8B9wOeBd4H/EpHnjDFbi/F5ilJp5Op56E90U5XSG0h5vUq/0alcCm+BhllPimWBnga0G2PeMsYcAX4OnFOkz1KUiiNVDdlCiVY612u2x1fxrGSKYoGGVk+KNQd6ArAz8Pxd4PTgDiJyJXClfTaySMMofyrJVaYXZp9UwpksfJncq9l+l7n+vvRvVOkUZQ60Xz0pFWKMKfxBRRYBf2GMucI9vxQ4zRjz12n2/wPwIbC34IMpLqOJ3pghmuOO4phBxz2YRHHMUJpxf9oY86lCH1REVmPPJ1eqgJ7A8weNMQ+6Y+akJ4NJsSzQd4G6wPMTgY50OxtjPiUiLcaYpiKNpyhEccwQzXFHccyg4x5MojhmiO64U2GMmV+Ew+akJ4NJseZA/wtoFJEJIjIEuAh4rkifpSiKopQvodWToligxpijInId8P8CHwN+bIz5bTE+S1EURSlfwqwnRSukYIx5Hng+h7c8WKyxFJEojhmiOe4ojhl03INJFMcM0R33oJGHngwKRQkiUhRFUZRyR2vhKoqiKEoelFxAw1qiKRUiEhORN0Rks4i0uG3VIvJrEWlzj58MwTh/LCJdItIa2JZynGL5V/f9vy4i00M05qUisst935tFZEHgtW+4Mb8pIn9RojHXichLIrJNRH4rIl9128P+Xacbd9i/7yoReUVEtrhx/4PbPkFENrnv+xcu0AQRGeqet7vX60M05kdEZEfgu57mtofiN6JkiTGmZAt2Qvj3wERgCLAFmFLKMfUz3hgwOmnbPwE3u/Wbge+FYJx/DkwHWvsbJ7AAeAEQYAawKURjXgrcmGLfKe63MhSY4H5DHyvBmMcB0936J4DfubGF/btON+6wf98CHOfWPw5sct/jE8BFbvsDwNVu/RrgAbd+EfCLEI35EeCCFPuH4jeiS3ZLqS3Q0JZoyoFzgJ+49Z8A55ZwLAAYY/4T2Je0Od04zwF+aiwbgVEiMm5wRuqTZszpOAf4uTGm1xizA2jH/pYGFWPMbmPMa279A2AbtmpK2L/rdONOR1i+b2OMOeieftwtBpgDrHTbk79v7++wEpgrIjJIwwUyjjkdofiNKNlRagFNVaIp0z9yqTHAf4jIq2JLEQKMNcbsBnthAsaUbHSZSTfOsP8NrnOurB8H3OOhG7NzD56CtTAi810njRtC/n2LyMdEZDPQBfwaaw3vN8Z49eOCY4uP273eDRw/uCPuO2ZjjPddf8d91/8sIkOTx+wo+W9ESU+pBTTV3WCYw4JnGWOmA2cB14rIn5d6QAUgzH+D+4E/BaYBu4G73PZQjVlEjgOeAq43xhzItGuKbWEad+i/b2PMR8aYadhqNKcBk1Pt5h5DMe7kMYvIVOAbwCTg/wKqgZvc7qEYs5IdpRbQ0JZoSoUxpsM9dgH/jv0H3uO5WNxjV+lGmJF04wzt38AYs8ddfP4I/AjfbRiaMYvIx7Ei9Jgx5mm3OfTfdapxR+H79jDG7AeasfOEo0TEy2kPji0+bvf6SLKfJig4gTHPd250Y4zpBR4mxN+1kp5SC2hoSzQlIyLHisgnvHXgfwGt2PF+xe32FeDZ0oywX9KN8zngL1303wyg23M/lpqkuZ/zsN832DFf5KIsJwCNwCslGJ8ADwHbjDF3B14K9XedbtwR+L4/JSKj3PowYB52/vYl4AK3W/L37f0dLgDWGmMG1ZpLM+btgRsswc7ZBr/rkv9GlCwpdRQTNursd9i5jFtLPZ4M45yIjUTcAvzWGyt2TuVFoM09VodgrD/DuuD+G3tHe3m6cWJdRve57/8NoClEY37Ujel17IVlXGD/W92Y3wTOKtGYz8S6114HNrtlQQS+63TjDvv3/VngN258rcDtbvtErKC3A08CQ932Kve83b0+MURjXuu+61ZgBX6kbih+I7pkt2glIkVRFEXJg1K7cBVFURQlkqiAKoqiKEoeqIAqiqIoSh6ogCqKoihKHqiAKoqiKEoeqIAqiqIoSh6ogCqKoihKHqiAKoqiKEoe/P97D0i6C4DhgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sample every 10th data point of the 'sic' variable\n",
    "dset = Dataset(outfile,'r')\n",
    "sicmap = dset.variables['sic'][::10,::10]\n",
    "\n",
    "# Make image with colour map jet (as example)\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.imshow(sicmap, origin='lower', cmap='jet') \n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('%')\n",
    "plt.title(dset.title)\n",
    "#plt.xlabel('X')\n",
    "#plt.ylabel('Y')\n",
    "image_file = 'tssic-mean.png' # input('Name of plot file (PNG): ')\n",
    "plt.savefig(image_file, bbox_inches=0)\n",
    "\n",
    "# Close data file\n",
    "dset.close()\n",
    "\n",
    "print('Graphics of mean sea ice concentration is stored in file', image_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}