diff --git a/README.md b/README.md index ba6bf09..82c2980 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,8 @@ Numerical models are widely used, but gaining expertise in how they work has oft - [JT Thielen](https://github.com/jthielen) - [Sam Gardner](https://github.com/wx4stg) - [Roger Riggin](https://github.com/riotrogerriot) + - [Justin Spotts](https://github.com/jrspotts) + - [Mathieu R](https://github.com/shenronUber) ### Contributors @@ -32,6 +34,10 @@ Numerical models are widely used, but gaining expertise in how they work has oft +Addition contributions to discussions and decisions for this notebook by: + +- [Rachel Smith](https://github.com/rachaellsmith) + ## Resources This cookbook would not be possible without the vast collection of academic texts and prior work in atmospheric modeling. The key resources used in building this notebook include: diff --git a/notebooks/config/grid_staggering.ipynb b/notebooks/config/grid_staggering.ipynb index 2441a9c..0eb2863 100644 --- a/notebooks/config/grid_staggering.ipynb +++ b/notebooks/config/grid_staggering.ipynb @@ -295,7 +295,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.8" + "version": "3.11.6" }, "nbdime-conflicts": { "local_diff": [ diff --git a/notebooks/images/basic_bubble.png b/notebooks/images/basic_bubble.png new file mode 100644 index 0000000..ffd01e1 Binary files /dev/null and b/notebooks/images/basic_bubble.png differ diff --git a/notebooks/images/c-grid_example.svg b/notebooks/images/c-grid_example.svg new file mode 100644 index 0000000..c58377b --- /dev/null +++ b/notebooks/images/c-grid_example.svg @@ -0,0 +1,1708 @@ + + + + + + + + 2024-06-14T12:56:13.472104 + image/svg+xml + + + Matplotlib v3.8.4, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/notebooks/intro/constants.py b/notebooks/intro/constants.py new file mode 100644 index 0000000..664c273 --- /dev/null +++ b/notebooks/intro/constants.py @@ -0,0 +1,7 @@ +#!/usr/bin/env python3 + +R_d = 8.314462618 / 28.96546e-3 # R / Md, (J/(mol*K)) / (kg/mol) = J/(kg*K) +c_p = 1.4 * R_d / (1.4 - 1) # J/(kg*K) +c_v = c_p - R_d +gravity = 9.80665 # m/(s^(2)) +p0 = 1.0e5 diff --git a/notebooks/intro/crash_course.ipynb b/notebooks/intro/crash_course.ipynb index 7bcc15a..d48e5c2 100644 --- a/notebooks/intro/crash_course.ipynb +++ b/notebooks/intro/crash_course.ipynb @@ -4,22 +4,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Let's start here! If you can directly link to an image relevant to your notebook, such as [canonical logos](https://github.com/numpy/numpy/blob/main/doc/source/_static/numpylogo.svg), do so here at the top of your notebook. You can do this with Markdown syntax,\n", - "\n", - "> `![](http://link.com/to/image.png \"image alt text\")`\n", - "\n", - "or edit this cell to see raw HTML `img` demonstration. This is preferred if you need to shrink your embedded image. **Either way be sure to include `alt` text for any embedded images to make your content more accessible.**\n", - "\n", - "\"Project" + "\"Colormesh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Crash Course Towards Your First Model\n", - "\n", - "Next, title your notebook appropriately with a top-level Markdown header, `#`. Do not use this level header anywhere else in the notebook. Our book build process will use this title in the navbar, table of contents, etc. Keep it short, keep it descriptive. Follow this with a `---` cell to visually distinguish the transition to the prerequisites section." + "# Crash Course Towards Your First Model" ] }, { @@ -34,13 +26,8 @@ "metadata": {}, "source": [ "## Overview\n", - "If you have an introductory paragraph, lead with it here! Keep it short and tied to your material, then be sure to continue into the required list of topics below,\n", "\n", - "1. This is a numbered list of the specific topics\n", - "1. These should map approximately to your main sections of content\n", - "1. Or each second-level, `##`, header in your notebook\n", - "1. Keep the size and scope of your notebook in check\n", - "1. And be sure to let the reader know up front the important concepts they'll be leaving with" + "The following notebook serves as a crash course in constructing a simple two-dimensional mesoscale atmospheric numerical model. To begin, we select a closed set of dynamical equations, in line with Klemp and Wilhelmson (1978; hereafter KW78). Necessary assumptions are stated to simplify the equations into a managble form to maximize both computational and educational applications. Pre-defined model configurations are presented (in an attempt) to ensure numerical stability. Equations are converted into finite-differences and then broken down into Python code to explicity demonstrate the construction of a numerical model. Additional resources will be provided in the future regarding varying model configurations and utilizing increasingly more realistic and complex model equations." ] }, { @@ -48,23 +35,14 @@ "metadata": {}, "source": [ "## Prerequisites\n", - "This section was inspired by [this template](https://github.com/alan-turing-institute/the-turing-way/blob/master/book/templates/chapter-template/chapter-landing-page.md) of the wonderful [The Turing Way](https://the-turing-way.netlify.app) Jupyter Book.\n", - "\n", - "Following your overview, tell your reader what concepts, packages, or other background information they'll **need** before learning your material. Tie this explicitly with links to other pages here in Foundations or to relevant external resources. Remove this body text, then populate the Markdown table, denoted in this cell with `|` vertical brackets, below, and fill out the information following. In this table, lay out prerequisite concepts by explicitly linking to other Foundations material or external resources, or describe generally helpful concepts.\n", - "\n", - "Label the importance of each concept explicitly as **helpful/necessary**.\n", "\n", "| Concepts | Importance | Notes |\n", "| --- | --- | --- |\n", - "| [Intro to Cartopy](https://foundations.projectpythia.org/core/cartopy/cartopy.html) | Necessary | |\n", - "| [Understanding of NetCDF](https://foundations.projectpythia.org/core/data-formats/netcdf-cf.html) | Helpful | Familiarity with metadata structure |\n", - "| Project management | Helpful | |\n", + "| [NumPy Basics](https://foundations.projectpythia.org/core/numpy/numpy-basics.html) | Necessary | |\n", + "| [Introduction to xarray](https://foundations.projectpythia.org/core/xarray/xarray-intro.html) | Helpful | Familiarity with coordinate-aware arrays |\n", + "| Dynamical Meteorology | Helpful | |\n", "\n", - "- **Time to learn**: estimate in minutes. For a rough idea, use 5 mins per subsection, 10 if longer; add these up for a total. Safer to round up and overestimate.\n", - "- **System requirements**:\n", - " - Populate with any system, version, or non-Python software requirements if necessary\n", - " - Otherwise use the concepts table above and the Imports section below to describe required packages as necessary\n", - " - If no extra requirements, remove the **System requirements** point altogether" + "- **Time to learn**: Estimated 30 to 60 minutes." ] }, { @@ -80,7 +58,7 @@ "source": [ "## Imports\n", "\n", - "Here we'll be using a basic set of Python libraries, along with Numba for fast numerical routines. A helpful constants file is also imported." + "Here we'll be using a basic set of Python libraries, along with Numba for fast numerical routines. A helpful constants file is also imported, along side the code that controls the running of the model. (For those curious, these files can be viewed [here](./constants) and [here](./constants).)" ] }, { @@ -94,111 +72,778 @@ "import xarray as xr\n", "import matplotlib.pyplot as plt\n", "\n", - "from constants import *" + "from constants import *\n", + "from driver import ModelDriver" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Introduction" + "## Basic Equations\n", + "\n", + "### Starting Equations (from WK78)\n", + "\n", + "#### Diagnostics\n", + " \n", + "DEQ 1. Equation of State (2.1) \n", + "$$\n", + "p = \\rho R_{d} T\n", + "$$ \n", + "* $p$: Pressure\n", + "* $\\rho$: Density\n", + "* $R_{d}$: Specific gas constant for dry air\n", + "* $T: Temperature\n", + "
\n", + "\n", + "DEQ 2. Exner Function (2.2)\n", + "$$\n", + "\\Pi = (\\frac{p}{p_{0}})^\\frac{R_{d}}{c_{p}}\n", + "$$\n", + "* $\\Pi$: Non-Dimensional Pressure\n", + "* $p$: Pressure\n", + "* $p_{0}$: Reference Pressure\n", + "* $R_{d}$: Specific gas constant for dry air\n", + "* $c_{p}$: Specific heat at constant pressure \n", + "
\n", + "\n", + "\n", + "#### Prognostics\n", + " \n", + "PEQ 1-2. Momentum Equation (Zonal & Vertical; 2.4)\n", + "$$\n", + "\\underset{1}{\\frac{du_{i}}{dt}}+ \n", + "\\underset{2}{c_{p} \\overline{\\theta_{v}} \\frac{\\partial \\pi}{\\partial x_{i}} }=\n", + "\\underset{3}{\\delta_{i_{3}}} \n", + "\\underset{4}{g}\\left[\n", + "\\underset{5}{\\frac{\\theta}{\\overline{\\theta}} - 1}+\n", + "\\underset{6}{0.61(q_{v}-\\overline{q_{v}})}-\n", + "\\underset{7}{q_{c}-q_{r}}\\right]-\n", + "\\underset{8}{\\epsilon_{ij_{3}}}fu_{i}+ \n", + "\\underset{9}{D_{u_i}}\n", + "$$\n", + "\n", + "1. Lagrangian of Wind\n", + "2. PGF Term\n", + "3. Kronecker Delta (i.e., the term that follows only appears form dimension 3, the vertical)\n", + "4. Gravity\n", + "5. Dry Buoyancy Contribution\n", + "6. Moist Buoyancy Contribution\n", + "7. Water Loading\n", + "8. Coriolis Term\n", + "9. Turblence Term\n", + "\n", + "Derived via Navier-Stokes equations, along with DEQ1-2, Hydrostatic Function, and Linearized Pressure Term. Tensor Notation is used for simplicity.\n", + "
\n", + " \n", + "PEQ 3. Prognostic Equations (2.5)\n", + "$$\n", + "\\underset{1}{\\frac{d\\phi}{dt}}=\n", + "\\underset{2}{M_{\\phi}}+ \n", + "\\underset{3}{D_{\\phi}}\n", + "$$\n", + "\n", + "1. Lagrangian of Prognositc Variable\n", + "2. Microphysical Term \n", + "3. Turbulence Term\n", + " \n", + "$\\phi$ is representative of either $\\theta, q_{v}, q_{r},$ or $q_{c}$\n", + "
\n", + " \n", + "PEQ 4. Pressure Equation (2.7a)\n", + "$$\n", + "\\underset{1}{\\frac{\\partial\\pi}{\\partial t}}+\n", + "\\underset{2}{\\frac{ \\overline{c}^2}{c_{p}{\\overline\\rho\\overline{\\theta_{v}^2}}}\n", + "{\\frac{\\partial}{\\partial x_{j}}(\\overline{\\rho}\\overline{\\theta_{v}}} u_{j})}=\n", + "\\underset{3}{f_{\\pi}}\n", + "$$\n", + "\n", + "1. Eularian of Pressure\n", + "2. Pressure Adjustment Term\n", + "3. Non-Relevant Terms for Sound & Gravity Waves (See KW78 2.7b)\n", + "\n", + "Derived using Compressible Continuity & Thermodynamic Equations; Tensor Notation Used for Simplicity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "...TODO..." + "### Assumptions and Simplification\n", + "\n", + "In this notebook we construct a two-dimensional, dry, and compressible atmospheric model that is broadly in line with KW78, though several assumptions and choice configurations were made to simplify the current model for computational and educational efficency.\n", + "
\n", + "1. We will only consider the zonal ($x$) and vertical ($z$) components.\n", + "2. Base-state variables are a function of $z$ only, denoted by $\\overline{\\phi}$.\n", + "3. Water-vapor is neglected (i.e, $q_{v}, q_{r}, q_{c} = 0$), so $T_{v} = T$ and/or $\\theta_{v} = \\theta$.\n", + "4. Coriolis force, microphysics, and Turbulence are also neglected(i.e., $f, M_{\\theta}, D_{\\theta} = 0$). \n", + "5. As in KW78, the $f_{\\pi}$ term in Pressure Equation (PEQ4-3) is neglected due to its negligible influences on convective-scale processes along with sound and gravity waves.\n", + "6. Sub-Grid Processes require parameterizations in order to achieve model closure. For example, sub-grid turbulence is first obtained using Reynolds Averaged prognostic variables (i.e., breaking up variables into mean and turbulence components), and then must be parameterized using additional assumptions (such as the flux-gradient approach). \n", + "7. The current test case is designed to have a calm and isentropic base-state (i.e., $\\frac{\\partial \\overline\\theta}{\\partial t}$ and $\\overline{U} = 0$)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Final Prognostic Equations\n", + "\n", + "The afformentioned assumptions allowed for the derivation of the following simplified equations that serve as the foundation for the numerical model.\n", + "
\n", + "\n", + "EQ1. Zonal Momentum Equation\n", + "$$\n", + "\\underset{1}{\\frac{du}{dt}}+\n", + "\\underset{2}{c_{p} \\bar{\\theta} \\frac{\\partial \\pi}{\\partial x}}\n", + "=0\n", + "$$\n", + "\n", + "1. Lagrangian Derivative of Zonal Wind\n", + "2. PGF Term\n", + "
\n", + "\n", + "EQ2. Vertical Momentum Equation\n", + "$$\n", + "\\underset{1}{\\frac{dw}{dt}}+\n", + "\\underset{2}{c_{p} \\bar{\\theta} \\frac{\\partial \\pi}{\\partial z}}= \n", + "\\underset{3}{g\\left[\\frac{\\theta}{\\bar{\\theta}} - 1\\right]}\n", + "$$\n", + "\n", + "1. Lagrangian Derivative of Vertical Wind\n", + "2. PGF Term\n", + "3. Dry Buoyancy Contribution\n", + "
\n", + "\n", + "EQ3. Prognostic Equations\n", + "$$\n", + "\\underset{1}{\\frac{d\\theta}{dt}}\n", + "=0\n", + "$$\n", + "\n", + "1. Lagrangian Derivative of Potential Temperature\n", + "
\n", + "\n", + "EQ4. Pressure Equation \n", + "$$\n", + "\\underset{1}{\\frac{\\partial\\pi}{\\partial t}}+ \n", + "\\underset{2}{\\frac{\\overline{c}^2}{c_{p}\\overline{\\rho}\\overline{\\theta}^2}[\n", + "\\frac{\\partial}{\\partial x}(\\overline{\\rho}\\overline{\\theta}u) +\n", + "\\frac{\\partial}{\\partial z}(\\overline{\\rho}\\overline{\\theta}w)]} \n", + "=0\n", + "$$\n", + "\n", + "1. Eularian Derivative of Pressure\n", + "2. Pressure Adjustment Terms" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model Configuration\n", + "\n", + "Summary: \n", + "- The current test case was configured using a 32 x 20 $km^2$ grid, including by 160 zonal grid points (nx) and 100 vertical grid points (nz), with both the horizontal and vertical grid-spacing ($\\Delta x, \\Delta z$) set to 200 m. Equations are integrated using a 0.1 s time-step ($\\Delta t$). \n", + "
\n", + "\n", + "Domain (32 km x 20 km):\n", + "- nx = 160\n", + "- nz = 100\n", + "- $dx, dz = 200 m$\n", + "
\n", + "\n", + "Grid Type (Staggered Grid: C) (In/Around Box-Good for Advection)\n", + "- Mass: Centered (i,k)\n", + "- Velocity: Edges \n", + " - u (i +/- 1/2, k)\n", + " - w (i, k +/- 1/2)\n", + "
\n", + "\n", + "TODO: Place Grid Visual Aid Here\n", + "
\n", + "\n", + "Boundary Conditions\n", + "- Free-Slip Lower\n", + "- Rigid Top \n", + "- Periodic Lateral \n", + "
\n", + "\n", + "Initial Conditions\n", + "- P_surf = 950 hPa, with remainder of atmosphere determined via hydrostatic-balance\n", + "- Theta = 300 K (isentropic)\n", + "- Calm (u, w) = 0\n", + "- CI: +3 K Warm Bubble $\\theta^{\\prime}$, with radius of 4 $km$ centered at $z = 2$ $km$. $p^{\\prime}$ (and thus $\\pi$) is adjusted accordingly.\n", + "
\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discritization\n", + "\n", + "### General Approach\n", + "\n", + "The following space/time discritization methods were used to convert the model equations into a code-ready format.\n", + "\n", + "Centered Spatial Differencing, on a C-grid:\n", + "$$\n", + "\\frac{\\partial \\phi}{\\partial x_i} \\approx \\frac{\\phi_{i+\\frac{1}{2}} - \\phi_{i-\\frac{1}{2}}}{\\Delta x_i}\n", + "$$ \n", + "\n", + "However, due to averaging, many of the formulations became equivalent to centered differencing on a non-staggered grid:\n", + "\n", + "$$\n", + "\\frac{\\partial \\phi}{\\partial x_i} \\approx \\frac{\\phi_{i+1} - \\phi_{i-1}}{2 \\Delta x_i}\n", + "$$ \n", + "\n", + "Leap-Frog Time Differencing:\n", + "$$\n", + "\\frac{\\partial \\phi}{\\partial t} \\approx \\frac{\\phi^{\\tau+1} - \\phi^{\\tau-1}}{2 \\Delta t}\n", + "$$\n", + "\n", + "### Equation-by-Equation\n", + "\n", + "#### u-momentum\n", + "\n", + "Given our indexing notation, this equation is centered on the $(i + \\frac{1}{2}, k)$ point.\n", + "\n", + "u advection Term:\n", + "$$\n", + "u\\frac{\\partial u}{\\partial x} \\cong\n", + "\\frac{1}{2\\Delta x}u_{i + \\frac{1}{2}, k}(u_{i + \\frac{3}{2}, k}-\n", + "u_{i - \\frac{1}{2}, k})\n", + "$$\n", + "\n", + "TODO: show full derivation, not just final equation\n", + "\n", + "To implement this in code, we must consider that the array indexes are whole numbers, and relative to the particular array. When we loop over the two spatial dimensions for $u$, this becomes:" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "# some subsection code\n", - "new = \"helpful information\"" + "@numba.njit()\n", + "def u_tendency_u_advection_term(u, dx):\n", + " term = np.zeros_like(u)\n", + " for k in range(1, u.shape[0] - 1):\n", + " for i in range(1, u.shape[1] - 1):\n", + " term[k, i] = u[k, i] * (u[k, i + 1] - u[k, i - 1]) / (2 * dx)\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Basic Equations\n", + "w advection term:\n", "\n", - "...TODO..." + "$$\n", + "w\\frac{\\partial u}{\\partial z} \\cong\n", + "\\frac{1}{4}( w_{i + 1, k + \\frac{1}{2}} + w_{i + 1, k - \\frac{1}{2}} +\n", + "w_{i, k + \\frac{1}{2}} + w_{i, k - \\frac{1}{2}})\n", + "\\frac{1}{2\\Delta z}(u_{i + \\frac{1}{2}, k + 1} - u_{i + \\frac{1}{2}, k - 1} )\n", + "$$\n", + "\n", + "TODO: show full derivation, not just final equation" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def u_tendency_w_advection_term(u, w, dz):\n", + " term = np.zeros_like(u)\n", + " for k in range(1, u.shape[0] - 1):\n", + " for i in range(1, u.shape[1] - 1):\n", + " term[k, i] = 0.25 * (\n", + " w[k + 1, i] + w[k + 1, i - 1] + w[k, i] + w[k, i - 1] \n", + " ) * (u[k + 1, i] - u[k - 1, i]) / (2 * dz)\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Starting Equations\n", + "PGF Term:\n", "\n", - "...TODO..." + "$$\n", + "c_{p}\\bar{\\theta} \\frac{\\partial \\pi}{\\partial x} \\cong \n", + "c_{p}\\bar{\\theta_{k}}\\frac{1}{\\Delta x}(\\pi_{i+1,k}-\\pi_{i,k} )\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def u_tendency_pgf_term(u, pi, theta_base, dx):\n", + " term = np.zeros_like(u)\n", + " for k in range(1, u.shape[0] - 1):\n", + " for i in range(1, u.shape[1] - 1):\n", + " term[k, i] = c_p * theta_base[k] * (pi[k, i] - pi[k, i - 1]) / dx\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Assumptions and Simplification\n", + "Now, we can combine all these RHS terms, accounting for the negation/subtraction present in the full equation" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def u_tendency(u, w, pi, theta_base, dx, dz):\n", + " return (\n", + " u_tendency_u_advection_term(u, dx)\n", + " + u_tendency_w_advection_term(u, w, dz)\n", + " + u_tendency_pgf_term(u, pi, theta_base, dx)\n", + " ) * -1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### w momentum\n", + "\n", + "Here, we proceed similarly\n", "\n", - "...TODO..." + "u advection term:\n", + "\n", + "$$\n", + "u\\frac{\\partial w}{\\partial x} = \n", + "\\frac{1}{4}(u_{i + \\frac{1}{2}, k + 1} +\n", + "u_{i + \\frac{1}{2}, k} + \n", + "u_{i - \\frac{1}{2}, k +1} +\n", + "u_{i - \\frac{1}{2}, k})\n", + "\\frac{1}{2\\Delta x}(w_{i + 1, k + \\frac{1}{2}} - \n", + "w_{i - 1, k + \\frac{1}{2}})\n", + "$$\n", + "\n", + "TODO: full derivation" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def w_tendency_u_advection_term(u, w, dx):\n", + " term = np.zeros_like(w)\n", + " for k in range(2, term.shape[0] - 2):\n", + " for i in range(1, term.shape[1] - 1):\n", + " term[k, i] = 0.25 * (\n", + " u[k, i] + u[k, i + 1] + u[k - 1, i] + u[k - 1, i + 1] \n", + " ) * (w[k, i + 1] - w[k, i - 1]) / (2 * dx)\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Final Prognostic Equations\n", + "w advection term:\n", + "\n", + "$$\n", + "w\\frac{\\partial w}{\\partial z}=\n", + "w_{i, k + \\frac{1}{2}}\n", + "\\frac{1}{2\\Delta z} (w_{i, k + \\frac{3}{2}} - w_{i-1, \\frac{1}{2}})\n", + "$$\n", "\n", - "...TODO..." + "TODO: full derivation" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def w_tendency_w_advection_term(w, dz):\n", + " term = np.zeros_like(w)\n", + " for k in range(2, term.shape[0] - 2):\n", + " for i in range(1, term.shape[1] - 1):\n", + " term[k, i] = w[k, i] * (w[k + 1, i] - w[k - 1, i]) / (2 * dz)\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Model Configuration\n", + "PGF term:\n", "\n", - "...TODO..." + "$$ \n", + "c_{p}\\overline{\\theta}\\frac{\\partial \\pi}{\\partial z} =\n", + "c_{p} \\frac{1}{2}( \\overline{\\theta}_{k+1} + \\overline{\\theta}_{k})\n", + "\\frac{1}{\\Delta z}(\\pi_{i, k+1} - \\pi_{i,k})\n", + "$$\n", + "\n", + "TODO: full derivation" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def w_tendency_pgf_term(w, pi, theta_base, dz):\n", + " term = np.zeros_like(w)\n", + " for k in range(2, term.shape[0] - 2):\n", + " for i in range(1, term.shape[1] - 1):\n", + " term[k, i] = c_p * 0.5 * (theta_base[k] + theta_base[k - 1]) * (pi[k - 1, i] - pi[k, i]) / dz\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Grid Type\n", + "Buoyancy Term:\n", + "\n", + "$$\n", + " g\\left[\\frac{\\theta^{\\prime}}{\\overline{\\theta}}\\right] =\n", + " g\\left[\\frac{\\theta_{i, k+1}^{\\prime} + \\theta_{i, k}^{\\prime}}{\\overline{\\theta}_{k+1}, \\overline{\\theta}_{k}}\\right]\n", + "$$\n", "\n", - "...TODO..." + "TODO: full derivation" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def w_tendency_buoyancy_term(w, theta_p, theta_base):\n", + " term = np.zeros_like(w)\n", + " for k in range(2, term.shape[0] - 2):\n", + " for i in range(1, term.shape[1] - 1):\n", + " term[k, i] = gravity * (theta_p[k, i] + theta_p[k - 1, i]) / (theta_base[k] + theta_base[k - 1]) \n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Initial/Boundary Conditions\n", + "Which, in combination, becomes:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def w_tendency(u, w, pi, theta_p, theta_base, dx, dz):\n", + " return (\n", + " w_tendency_u_advection_term(u, w, dx) * -1.0\n", + " - w_tendency_w_advection_term(w, dz)\n", + " - w_tendency_pgf_term(w, pi, theta_base, dz)\n", + " + w_tendency_buoyancy_term(w, theta_p, theta_base)\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Theta tendency equation\n", + "\n", + "u advection term:\n", + "\n", + "$$\n", + " u\\frac{\\partial \\theta^{\\prime}}{\\partial x} = \n", + " \\frac{1}{2}(u_{i+\\frac{1}{2}, k} + u_{i-\\frac{1}{2}})\n", + " \\frac{1}{2\\Delta x}(\\theta_{i+1,k}^{\\prime}-\\theta_{i-1,k}^{\\prime})\n", + " $$" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def theta_p_tendency_u_advection_term(u, theta_p, dx):\n", + " term = np.zeros_like(theta_p)\n", + " for k in range(1, theta_p.shape[0] - 1):\n", + " for i in range(1, theta_p.shape[1] - 1):\n", + " term[k, i] = (\n", + " (u[k, i + 1] + u[k, i]) / 2\n", + " * (theta_p[k, i + 1] - theta_p[k, i - 1]) / (2 * dx)\n", + " )\n", + " return term" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "w advection of theta perturbation:\n", "\n", - "...TODO..." + "$$\n", + " w\\frac{\\partial \\theta^{\\prime}}{\\partial z} = \n", + " \\frac{1}{2}(w_{i+\\frac{1}{2}, k} + w_{i-\\frac{1}{2}})\n", + " \\frac{1}{2\\Delta z}(\\theta_{i,k+1}^{\\prime}-\\theta_{i,k-1}^{\\prime})\n", + " $$" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def theta_p_tendency_w_advection_of_perturbation_term(w, theta_p, dz):\n", + " term = np.zeros_like(theta_p)\n", + " for k in range(1, theta_p.shape[0] - 1):\n", + " for i in range(1, theta_p.shape[1] - 1):\n", + " term[k, i] = (\n", + " (w[k + 1, i] + w[k, i]) / 2\n", + " * (theta_p[k + 1, i] - theta_p[k - 1, i]) / (2 * dz)\n", + " )\n", + " return term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Discritization\n", + "w advection of theta base state:\n", + "\n", + "$$\n", + " w\\frac{\\partial \\overline{\\theta}}{\\partial z} = \n", + " \\frac{1}{2}(w_{i, k+\\frac{1}{2}} + w_{i,k\\frac{1}{2}})\n", + " \\frac{1}{2\\Delta z}(\\overline{\\theta}_{k+1}-\\overline\\theta_{k-1})\n", + " $$" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def theta_p_tendency_w_advection_of_base_term(w, theta_p, theta_base, dz):\n", + " term = np.zeros_like(theta_p)\n", + " for k in range(1, theta_p.shape[0] - 1):\n", + " for i in range(1, theta_p.shape[1] - 1):\n", + " term[k, i] = (\n", + " (w[k + 1, i] + w[k, i]) / 2\n", + " * (theta_base[k + 1] - theta_base[k - 1]) / (2 * dz)\n", + " )\n", + " return term" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Combining, becomes:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def theta_p_tendency(u, w, theta_p, theta_base, dx, dz):\n", + " return (\n", + " theta_p_tendency_u_advection_term(u, theta_p, dx)\n", + " + theta_p_tendency_w_advection_of_perturbation_term(w, theta_p, dz)\n", + " + theta_p_tendency_w_advection_of_base_term(w, theta_p, theta_base, dz)\n", + " ) * -1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Non-dimensional Pressure Tendency\n", + "\n", + "Based on the treatment by KW78, this equation expressed with a leading factor (which has $k$-index dependence) and two interior terms (which have finite differences). We express these together as:\n", + "\n", + "\\begin{align*}\n", + "-\\frac{\\bar{c}^2}{c_p \\bar{\\rho}_k \\bar{\\theta}_k^2}\\bigg[ &\\bar{\\rho}_k \\bar{\\theta}_k \\frac{1}{\\Delta x} \\left(u_{i + \\frac{1}{2}, k} - u_{i - \\frac{1}{2}, k}\\right) \\\\\n", + "\\quad\\quad &+ \\frac{1}{2}\\left(w_{i, k+\\frac{1}{2}} + w_{i, k-\\frac{1}{2}}\\right) \\frac{1}{2\\Delta z} \\left(\\bar{\\rho}_{k+1} \\bar{\\theta}_{k+1} - \\bar{\\rho}_{k-1} \\bar{\\theta}_{k-1}\\right) \\\\\n", + "\\quad\\quad &+ \\bar{\\rho}_k \\bar{\\theta}_k \\frac{1}{\\Delta z} \\left(w_{i , k+ \\frac{1}{2}} - w_{i, k - \\frac{1}{2}}\\right)\\bigg]\n", + "\\end{align*}" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "@numba.njit()\n", + "def pi_tendency(u, w, pi, theta_base, rho_base, c_s_sqr, dx, dz):\n", + " term = np.zeros_like(pi)\n", + " for k in range(1, pi.shape[0] - 1):\n", + " for i in range(1, pi.shape[1] - 1):\n", + " term[k, i] = (\n", + " -1 * (c_s_sqr / (c_p * rho_base[k] * theta_base[k]**2))\n", + " * (\n", + " (rho_base[k] * theta_base[k] * (u[k, i + 1] - u[k, i]) / dx)\n", + " + (\n", + " (w[k + 1, i] + w[k, i]) / 2\n", + " * (rho_base[k + 1] * theta_base[k + 1] - rho_base[k - 1] * theta_base[k - 1]) / (2 * dz)\n", + " )\n", + " + (rho_base[k] * theta_base[k] * (w[k + 1, i] - w[k, i]) / dz)\n", + " )\n", + " )\n", + " return term" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Running the Model\n", + "\n", + "We are now ready to set up and run the model!\n", + "\n", + "### Setup\n", + "\n", + "To start, we initialize the model driver class (see the [`driver.py`](./driver) if interested in the details) using our aforementioned model options, and plugging in our tendency equations written above). We then add in the initial base state and perturbations, so that the model has something to simulate. Finally, we export the initial state to inspect later:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up model\n", + "model = ModelDriver(\n", + " nx=160, nz=100, dx=200., dz=200., dt=0.01, c_s_sqr=50.0**2,\n", + " u_tendency=u_tendency, w_tendency=w_tendency, theta_p_tendency=theta_p_tendency, pi_tendency=pi_tendency\n", + ")\n", + "model.initialize_isentropic_base_state(300., 9.5e4)\n", + "model.initialize_warm_bubble(3.0, 4.0e3, 4.0e3, 2.0e3)\n", + "\n", + "# Start saving results\n", + "results = []\n", + "results.append(model.current_state())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can quickly preview what this initial warm bubble looks like using xarray:" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 51, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAHFCAYAAAD40125AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABaGklEQVR4nO3de1zUVf4/8NdnBmZAhFFAGEhE7UIqVIatoJuXLoCJrtWmpl9+Wka1pmbqttlN7VtaWWZb282veUlLv6VWrsaqpbaFl0T9LqaZliYqiCnMCCIDM+f3xzCf+YyDwDAzMDO8no/H5+FwPudz5pwp4M25SkIIASIiIiICAKhauwJEREREvoTBEREREZECgyMiIiIiBQZHRERERAoMjoiIiIgUGBwRERERKTA4IiIiIlJgcERERESkwOCIiIiISIHBEZGfy8/Px+zZs1FeXt4q719RUYGpU6ciPj4eISEhuOmmm7Bq1aomP19aWorx48cjOjoa7dq1Q3p6Or7++usmPdu1a1dIkgRJkjBp0iSX6v3555/Lz0qShD179rj0PBEFLgZHRH4uPz8fc+bMabXg6J577sGyZcswa9YsfPXVV7jllltw//334+OPP2702erqatx+++34+uuv8eabb+KLL75AbGwssrKysH379ia9/1133YUdO3ZgxowZLtV74MCB2LFjB5599lmXniOiwBfU2hUgIv+1ceNGbN68GR9//DHuv/9+AMDgwYPx22+/4a9//StGjRoFtVp9xecXL16MAwcOID8/H+np6fLzN954I5588kns2rWr0Tp06tQJaWlpLte9Y8eOSEtLw08//eTys0QU2NhzROTHZs+ejb/+9a8AgG7duslDRNu2bWuR91+3bh3at2+P++67zyH9gQcewOnTpxsNbtatW4ekpCQ5MAKAoKAg/Nd//Rd2796NU6dONateFosFL774IpKSkhAaGooOHTrghhtuwJtvvtms8oiobWHPEZEfe+ihh3D+/Hm89dZbWLt2LeLi4gAAPXv2vOIzQgiYzeYmlR8U1PCPiAMHDqBHjx5O+W644Qb5fr9+/Rp8/tZbb3VKtz3/448/4qqrrmpSXZVeffVVzJ49G88++ywGDBiAmpoa/PTTT6029EhE/oXBEZEf69y5M7p06QIA6N27N7p27droM8uWLcMDDzzQpPKFEA3eP3fuHLp37+6UHhkZKd9v7Hlb3uY8fyXff/89UlJSMHv2bDktMzOzWWURUdvD4IiojRk2bBh++OEHj5UnSVKz7nnq+fr84Q9/wIYNGzBx4kT86U9/Qnp6OiIiIppVFhG1PQyOiNqYyMhI6HQ6j5QVFRVVb+/O+fPn5ffy5vNXMnPmTISFhWHFihV47733oFarMWDAALzyyivo06dPs8okoraDE7KJ2phly5YhODi4SVdjUlJScOjQIdTW1jqkFxYWAgCSk5Mbfd6WtznPX0lQUBCmTZuGvXv34vz58/jkk09QVFSEzMxMXLx4sVllElHbwZ4jIj+n1WoBAFVVVU3K78lhtbvvvhuLFi3CmjVrMGrUKDl92bJliI+PR9++fRt9fuLEidi1a5ect7a2FitWrEDfvn0RHx/vdh07dOiAP//5zzh16hSmTp2K48ePNzhhnYiIwRGRn0tJSQEAvPnmmxg3bhyCg4ORlJSE8PDwevNHRUUhKirKI+89ZMgQ3HnnnfjLX/4Co9GIa665Bp988gny8vKwYsUKhz2OJkyYgGXLluGXX35BYmIiAODBBx/EP/7xD9x33314+eWXERMTg3feeQeHDx/Gli1bml2vYcOGITk5GX369EGnTp3w22+/YeHChUhMTMS1117rdruJKLAxOCLyc4MGDcLMmTOxbNkyLFq0CBaLBVu3bsWgQYNa5P3Xrl2LZ555Bs8//zzOnz+P66+/Hp988glGjx7tkM9sNsNsNjusgNNqtfj666/x5JNPYvLkybh48SJuuukmfPXVVxg4cGCz6zR48GCsWbMG//M//wOj0Qi9Xo8777wTzz33XJOGC4mobZNEY2t1iYh8VNeuXTFw4EAsXrwYKpUKKlXTp1Ha9ntavnw5JkyYgB9++IGTtYkIACdkE5GfW758OYKDgzFlyhSXnvviiy8QHByMCRMmeKlmROSv2HNERH6rsLAQ1dXVAICYmBh5Q8ymKC8vx9GjR+Wve/bsiXbt2nm8jkTkfxgcERERESlwWI2IiIg86t1338UNN9yAiIgIREREID09HV999VVrV6vJ2HNEREREHrV+/Xqo1Wpcc801AKx7n82fPx/79u1Dr169Wrl2jWNwRERERF4XGRmJ+fPn+8UiCO5z5EEWiwWnT59GeHh4sw/MJCKiwCeEwIULFxAfH+/SFhSuunTpEkwmk0fKEkI4/W7TarXyLv1XYjab8emnn6KyshLp6ekeqYu3MTjyoNOnTyMhIaG1q0FERH6iqKgInTt39krZly5dQlRoe1yE2SPltW/fHhUVFQ5ps2bNwuzZs+vNX1hYiPT0dFy6dAnt27fHunXr/OboHgZHHmQ7ruHokSNXPLqBiIjowoULuObaa736u8JkMuEizBiLq6Bxc/2VCRasrDiFoqIiREREyOkN9RolJSVh//79KC8vx5o1azBu3Dhs377dLwIkBkceZOtuDA8Pd/ifh4iIqD4tMQUjFCpoJPeCI3Xd7GTb6rOm0Gg08oTsPn364IcffsCbb76J999/3626tAQGR0RERAFMLUlQuxmEqSEBbi7fEkLIm7b6OgZHREREAUwlAWo3O6hUgEvB0dNPP40hQ4YgISEBFy5cwKpVq7Bt2zbk5eW5V5EWwuCIiIiIPOrMmTPIyclBcXExdDodbrjhBuTl5eHOO+9s7ao1CYMjIiKiAOaxYTUXLF682K33a20MjoiIiAKY2gPDamrPVMVv8Gw1IiIiIgX2HBEREQWw1hhW83cMjoiIiAIYh9Vcx2E1IiIiIgX2HBEREQUwDqu5jsERERFRAJPg/jBR2wqNOKxGRERE5IA9R0RERAGMw2quY3BEREQUwLhazXUMjoiIiAKYNThyt+eobeGcIyIiIiIF9hwREREFMA6ruY7BERERUQDjhGzXcViNiIiISIE9R0RERAFM5YFhtbbWk8LgiIiIKIBxWM11bS0YJCIiImoQe46IiIgCGFeruY7BERERUQBjcOQ6DqsRERERKbDniIiIKIBxQrbrGBwREREFMDU8MKwmPFIVv8HgiIiIKICpPNBzpHLzeX/DOUdERERECuw5IiIiCmAeWa3WtjqOGBwREREFMo9MyOawGhEREVHbxZ4jIiKiAMZhNdcxOCIiIgpgHFZzHYfViIiIiBTYc0RERBTAVJLk9j5FbW2fIwZHREREAUxSS5BU7gU3UhsLjlp1WO3bb7/FsGHDEB8fD0mS8PnnnzvclySp3mv+/PlynkGDBjndHz16tEM5ZWVlyMnJgU6ng06nQ05ODsrLyx3ynDhxAsOGDUNYWBiio6MxZcoUmEwmbzWdiIiIfFSr9hxVVlbixhtvxAMPPIB7773X6X5xcbHD11999RUmTJjglDc3NxcvvPCC/HVoaKjD/TFjxuDkyZPIy8sDADz88MPIycnB+vXrAQBmsxlDhw5Fp06d8N133+HcuXMYN24chBB46623PNJWIiKi1qBSS1C52XPEYbUWNGTIEAwZMuSK9/V6vcPXX3zxBQYPHozu3bs7pLdr184pr82hQ4eQl5eHnTt3om/fvgCARYsWIT09HYcPH0ZSUhI2bdqEgwcPoqioCPHx8QCA119/HePHj8dLL72EiIgId5pJRETUetQqSCo3B4qktnXyrN+sVjtz5gw2bNiACRMmON1buXIloqOj0atXL8yYMQMXLlyQ7+3YsQM6nU4OjAAgLS0NOp0O+fn5cp7k5GQ5MAKAzMxMVFdXo6Cg4Ip1qq6uhtFodLiIiIh8iaSSrPOO3Lnc7HnyN34zIXvZsmUIDw/HPffc45A+duxYdOvWDXq9HgcOHMDMmTPxf//3f9i8eTMAoKSkBDExMU7lxcTEoKSkRM4TGxvrcL9jx47QaDRynvrMmzcPc+bMcbdpRERE5EP8Jjj68MMPMXbsWISEhDik5+bmyq+Tk5Nx7bXXok+fPti7dy9uvvlmAPXPshdCOKQ3Jc/lZs6ciWnTpslfG41GJCQkNL1RREREXqZSS1C5ucW1Cuw58jn//ve/cfjwYaxevbrRvDfffDOCg4Nx5MgR3HzzzdDr9Thz5oxTvrNnz8q9RXq9Hrt27XK4X1ZWhpqaGqceJSWtVgutVutia4iIiFqOpHJ/zpEkOOfI5yxevBipqam48cYbG837448/oqamBnFxcQCA9PR0GAwG7N69W86za9cuGAwG9OvXT85z4MABh9VxmzZtglarRWpqqodbQ0RERL6sVXuOKioqcPToUfnrY8eOYf/+/YiMjESXLl0AWIeqPv30U7z++utOz//yyy9YuXIl7rrrLkRHR+PgwYOYPn06evfujf79+wMAevTogaysLOTm5uL9998HYF3Kn52djaSkJABARkYGevbsiZycHMyfPx/nz5/HjBkzkJuby5VqRETk1zis5rpW7Tnas2cPevfujd69ewMApk2bht69e+P555+X86xatQpCCNx///1Oz2s0Gnz99dfIzMxEUlISpkyZgoyMDGzZsgVqtVrOt3LlSqSkpCAjIwMZGRm44YYb8NFHH8n31Wo1NmzYgJCQEPTv3x8jR47EiBEj8Nprr3mx9URERN7n9kq1uqstkYRoYwOJXmQ0GqHT6XCmpIQ9TkREdEVGoxGxej0MBoPXfl/YfietT0lFmKLDoDkqzWYMKyzwan19iV9MyCYiIqLmsfb8uDkhGxYP1cY/MDgiIiIKYJxz5Dq/WK1GRERE1FLYc0RERBTAJMn94z8kS9vqOWJwREREFMBUahVUbs45Uom2NdDE4IiIiCiAeWIpviTaVs9R2woFiYiIyKvmzZuHW265BeHh4YiJicGIESNw+PDh1q6WSxgcERERBbCW3gRy+/bteOyxx7Bz505s3rwZtbW1yMjIQGVlpRdb6VkcViMiIgpgLT3nKC8vz+HrJUuWICYmBgUFBRgwYIBb9Wgp7DkiIiIirzEYDACAyMjIVq5J07HniIiIKJB54my0ugnZRqPRIVmr1UKr1V75MSEwbdo0/PGPf0RycrJ7dWhB7DkiIiIKYCpJgkrl5iVZg6OEhATodDr5mjdvXoPvPWnSJPznP//BJ5980hJN9Rj2HBEREVGTFBUVORw821Cv0eTJk/Hll1/i22+/RefOnVuieh7D4IiIiCiASWqV+wfPWqzPR0REOARH9RFCYPLkyVi3bh22bduGbt26ufXerYHBERERUQDzyMGzLhwf8thjj+Hjjz/GF198gfDwcJSUlAAAdDodQkND3apHS+GcIyIiIvKYd999FwaDAYMGDUJcXJx8rV69urWr1mTsOSIiIgpgHjk+xIWeIyGEW+/lCxgcERERBTBPzjlqKxgcERERBTCVGh6Yc+ShyviJthUKEhERETWCPUdEREQBTFJJkFRuzjly83l/w+CIiIgogKlUHjh41ty2BpraVmuJiIiIGsGeIyIiogDmkaX87h5c62cYHBEREQUwjyzld/N5f9O2WktERETUCPYcERERBTBJpYKkcrPnyM3n/Q2DIyIiogCmUntgtRqH1YiIiIjaLvYcERERBTIPTMhGG+s5YnBEREQUwCSVB1arcc4RERERBQpOyHZd22otERERUSPYc0RERBTArJtAqt0sw+yh2vgHBkdEREQBjDtku65VW/vtt99i2LBhiI+PhyRJ+Pzzzx3ujx8/HpIkOVxpaWkOeaqrqzF58mRER0cjLCwMw4cPx8mTJx3ylJWVIScnBzqdDjqdDjk5OSgvL3fIc+LECQwbNgxhYWGIjo7GlClTYDKZvNFsIiIi8mGtGhxVVlbixhtvxNtvv33FPFlZWSguLpavjRs3OtyfOnUq1q1bh1WrVuG7775DRUUFsrOzYTbbuwDHjBmD/fv3Iy8vD3l5edi/fz9ycnLk+2azGUOHDkVlZSW+++47rFq1CmvWrMH06dM932giIqIWpFKpPHK1Ja06rDZkyBAMGTKkwTxarRZ6vb7eewaDAYsXL8ZHH32EO+64AwCwYsUKJCQkYMuWLcjMzMShQ4eQl5eHnTt3om/fvgCARYsWIT09HYcPH0ZSUhI2bdqEgwcPoqioCPHx8QCA119/HePHj8dLL72EiIgID7aaiIio5XBYzXU+39pt27YhJiYG1113HXJzc1FaWirfKygoQE1NDTIyMuS0+Ph4JCcnIz8/HwCwY8cO6HQ6OTACgLS0NOh0Ooc8ycnJcmAEAJmZmaiurkZBQYG3m0hEREQ+xKcnZA8ZMgT33XcfEhMTcezYMTz33HO47bbbUFBQAK1Wi5KSEmg0GnTs2NHhudjYWJSUlAAASkpKEBMT41R2TEyMQ57Y2FiH+x07doRGo5Hz1Ke6uhrV1dXy10ajsdltJSIi8gb2HLnOp4OjUaNGya+Tk5PRp08fJCYmYsOGDbjnnnuu+JwQApIkyV8rX7uT53Lz5s3DnDlzGm0HERFRa5EkD2wCKbWt4MivWhsXF4fExEQcOXIEAKDX62EymVBWVuaQr7S0VO4J0uv1OHPmjFNZZ8+edchzeQ9RWVkZampqnHqUlGbOnAmDwSBfRUVFbrWPiIiIWp9fBUfnzp1DUVER4uLiAACpqakIDg7G5s2b5TzFxcU4cOAA+vXrBwBIT0+HwWDA7t275Ty7du2CwWBwyHPgwAEUFxfLeTZt2gStVovU1NQr1ker1SIiIsLhIiIi8iW2YTV3r7akVYfVKioqcPToUfnrY8eOYf/+/YiMjERkZCRmz56Ne++9F3FxcTh+/DiefvppREdH4+677wYA6HQ6TJgwAdOnT0dUVBQiIyMxY8YMpKSkyKvXevTogaysLOTm5uL9998HADz88MPIzs5GUlISACAjIwM9e/ZETk4O5s+fj/Pnz2PGjBnIzc1lwENERH6Nc45c16rB0Z49ezB48GD562nTpgEAxo0bh3fffReFhYVYvnw5ysvLERcXh8GDB2P16tUIDw+Xn3njjTcQFBSEkSNHoqqqCrfffjuWLl0KtWKr9JUrV2LKlCnyqrbhw4c77K2kVquxYcMGTJw4Ef3790doaCjGjBmD1157zdsfARERkVep1Cqo3Axu3H3e30hCCNHalQgURqMROp0OZ0pK2ONERERXZDQaEavXw2AweO33he130pE5DyM8RONWWRcumXDtrA+8Wl9f4tOr1YiIiMg9kkpyf7Wa6sortwMRgyMiIqIAxjlHrmtbrSUiIiJqBHuOiIiIAhh7jlzH4IiIiCiAcYds17Wt1hIRERE1gj1HREREAUxSq6FS7P3X3DLaEgZHREREAYxzjlzXtlpLRERE1Aj2HBEREQUw9hy5jsERERFRAJNUHlit5ubz/obBERERUQBjz5Hr2lZriYiIiBrBniMiIqIAJqkk93uOePAsERERBQrOOXJd22otERERUSPYc0RERBTAJJUaksrNHbLdfN7fMDgiIiIKZCq19XK3jDaEw2pERERECuw5IiIiCmQqlfVyt4w2hMERERFRAJPUakhqN+ccufm8v2lboSARERFRIxgcERERBTLbhGx3Lxd8++23GDZsGOLj4yFJEj7//HPvtM1LGBwREREFMpXKA8GRa+FCZWUlbrzxRrz99tteapR3cc4RERFRAGuNHbKHDBmCIUOGuPWerYnBERERETWJ0Wh0+Fqr1UKr1bZSbbyHw2pERESBTPLAfCPJOucoISEBOp1OvubNm9fKjfMO9hwREREFMg/ukF1UVISIiAg5ORB7jQAGR0RERNREERERDsFRoGJwREREFMBaY0K2v2NwREREFMha4eDZiooKHD16VP762LFj2L9/PyIjI9GlSxf36tICGBwRERGRR+3ZsweDBw+Wv542bRoAYNy4cVi6dGkr1arpGBwREREFMtsmkO6W4YJBgwZBCOHee7YiBkdEREQBjAfPuq5tzbAiIiIiagR7joiIiAKZSuXysFi9ZbQhrdrahk7trampwd/+9jekpKQgLCwM8fHx+H//7//h9OnTDmUMGjQIkiQ5XKNHj3bIU1ZWhpycHHlHz5ycHJSXlzvkOXHiBIYNG4awsDBER0djypQpMJlM3mo6ERFRy3D70FkPrHbzM60aHDV0au/Fixexd+9ePPfcc9i7dy/Wrl2Ln3/+GcOHD3fKm5ubi+LiYvl6//33He6PGTMG+/fvR15eHvLy8rB//37k5OTI981mM4YOHYrKykp89913WLVqFdasWYPp06d7vtFEREQtSFKpPXK1Ja06rNbQqb06nQ6bN292SHvrrbfwhz/8ASdOnHDYJ6Fdu3bQ6/X1lnPo0CHk5eVh586d6Nu3LwBg0aJFSE9Px+HDh5GUlIRNmzbh4MGDKCoqQnx8PADg9ddfx/jx4/HSSy+1id1AiYiIyMqvBhENBgMkSUKHDh0c0leuXIno6Gj06tULM2bMwIULF+R7O3bsgE6nkwMjAEhLS4NOp0N+fr6cJzk5WQ6MACAzMxPV1dUoKCi4Yn2qq6thNBodLiIiIp8iqezzjpp7SX4VLrjNbyZkX7p0CU899RTGjBnj0JMzduxYdOvWDXq9HgcOHMDMmTPxf//3f3KvU0lJCWJiYpzKi4mJQUlJiZwnNjbW4X7Hjh2h0WjkPPWZN28e5syZ44nmEREReYUnhsU4rOaDampqMHr0aFgsFrzzzjsO93Jzc+XXycnJuPbaa9GnTx/s3bsXN998MwBAkiSnMoUQDulNyXO5mTNnyrt+AoDRaERCQkLTG0ZEREQ+x+f7yWpqajBy5EgcO3YMmzdvbnT+z80334zg4GAcOXIEAKDX63HmzBmnfGfPnpV7i/R6vVMPUVlZGWpqapx6lJS0Wq18QnFbOamYiIj8jG2HbLcunw8XPMqnW2sLjI4cOYItW7YgKiqq0Wd+/PFH1NTUIC4uDgCQnp4Og8GA3bt3y3l27doFg8GAfv36yXkOHDiA4uJiOc+mTZug1WqRmprq4VYRERG1IHfnG3linyQ/06rDag2d2hsfH48///nP2Lt3L/75z3/CbDbLvTuRkZHQaDT45ZdfsHLlStx1112Ijo7GwYMHMX36dPTu3Rv9+/cHAPTo0QNZWVnIzc2Vl/g//PDDyM7ORlJSEgAgIyMDPXv2RE5ODubPn4/z589jxowZyM3NZW8QERFRG9OqwVFDp/bOnj0bX375JQDgpptucnhu69atGDRoEDQaDb7++mu8+eabqKioQEJCAoYOHYpZs2ZBrTgHZuXKlZgyZQoyMjIAAMOHD3fYW0mtVmPDhg2YOHEi+vfvj9DQUIwZMwavvfaat5pORETUIni2musk4c/H5voYo9EInU6HMyUl7HEiIqIrMhqNiNXrYTAYvPb7wvY76fyufyKifZh7ZVVUIrJvtlfr60v8YrUaERERNZMnjv9oY0v529YMKyIiIqJGsOeIiIgogEkqFSQ3V5u5+7y/YXBEREQUyCQPDKtJbWtYjcERERER+ZXmnGXqykRyBkdERESBTJLcPzi2gaO0WkOHDh0aPN7rcpIk4eeff0b37t2blJ/BERERUSCTVB4IjnxvztFnn32GyMjIRvMJIXDXXXe5VDaDIyIiIvIriYmJGDBgQJOOFQOA7t27Izg4uMnlMzgiIiIKYEJSQbjZ8+Pu85527Ngxl/IfOHDApfy+1VoiIiLyLNuwmruXj2lKwPPyyy83q2zfay0RERFRIzIzM3H8+PEr3n/llVcwa9asZpXN4IiIiCiQSZJnLh9z66234s4770RpaanTvfnz5+O5557DihUrmlU2gyMiIqJAplJ55vIxK1aswDXXXIOMjAwYDAY5/fXXX8fTTz+N5cuX47777mtW2b7XWiIiIvIY24Rsdy9fExQUhLVr16J9+/bIzs7GpUuXsHDhQjz11FNYtmwZRo8e3fyyPVhPIiIiohYTGhqKDRs2YODAgUhNTcXPP/+MJUuWYMyYMW6Vy+CIiIgokAXoJpBffvml/Povf/kLHn/8cdx9992IiIhwuDd8+HCXy2ZwREREFMgCNDgaMWKEU9pnn32Gzz77TP5akiSYzWaXy2ZwRERERH7HYrF4rWwGR0RERIEsQHuOvKlttZaIiKiNEZLkgdVqvrXP0Zdffomampom59+4cSOqqqqanJ/BEREREfmVu+++G+Xl5U3OP3r0aBQXFzc5P4fViIiIAlkADqsJITB+/Hhotdom5b906ZJL5TM4IiIiCmSeOP7Dx4bVxo0b51L+sWPHIiIiosn5GRwRERGRX1myZIlXy29ScGQ0Gl0u2JUIjYiIiLwkAIfVvK1JwVGHDh0gudClJkkSfv75Z3Tv3r3ZFSMiIiL3eeJsNF88W82bmjys9tlnnyEyMrLRfEII3HXXXW5VioiIiDxEUgEq9hy5oknBUWJiIgYMGICoqKgmFdq9e3cEBwe7VTEiIiKi1tCk4OjYsWMuFXrgwIFmVYaIiIg8jHOOXNa2WktERNTW2IIjdy8f9vXXXyM7OxtXX301rrnmGmRnZ2PLli3NLq9ZS/l3796Nbdu2obS01OngtwULFjS7MkRERESuePvtt/HEE0/gz3/+Mx5//HEAwM6dO3HXXXdhwYIFmDRpkstluhwczZ07F88++yySkpIQGxvrsIrNlRVtRERE1AICfFht3rx5eOONNxyCoClTpqB///546aWXWiY4evPNN/Hhhx9i/PjxLr8ZERERtSzbwbPuluGrjEYjsrKynNIzMjLwt7/9rVlluvxpqVQq9O/fv1lvRkRERORJw4cPx7p165zSv/jiCwwbNqxZZbrcc/TEE0/gH//4BxYuXNisNyQiIqIWFODDaj169MBLL72Ebdu2IT09HYB1ztH333+P6dOn4+9//7ucd8qUKU0qUxJCCFcqYbFYMHToUPz888/o2bOn035Ga9eubXJZ3377LebPn4+CggIUFxdj3bp1GDFihHxfCIE5c+bggw8+QFlZGfr27Yt//OMf6NWrl5ynuroaM2bMwCeffIKqqircfvvteOedd9C5c2c5T1lZGaZMmYIvv/wSgDXKfOutt9ChQwc5z4kTJ/DYY4/hm2++QWhoKMaMGYPXXnsNGo2mye0xGo3Q6XQ4U1LC41OIiOiKjEYjYvV6GAwGr/2+sP1OKj11wu33MBqNiLmqi1fr21zdunVrUj5JkvDrr782Ka/LPUeTJ0/G1q1bMXjwYERFRbk1CbuyshI33ngjHnjgAdx7771O91999VUsWLAAS5cuxXXXXYcXX3wRd955Jw4fPozw8HAAwNSpU7F+/XqsWrUKUVFRmD59OrKzs1FQUAC1Wg0AGDNmDE6ePIm8vDwAwMMPP4ycnBysX78eAGA2mzF06FB06tQJ3333Hc6dO4dx48ZBCIG33nqr2e0jIiIi73J1L8amcLnnKDw8HKtWrcLQoUM9WxFJcug5EkIgPj4eU6dOlSdUVVdXIzY2Fq+88goeeeQRGAwGdOrUCR999BFGjRoFADh9+jQSEhKwceNGZGZm4tChQ+jZsyd27tyJvn37ArB2t6Wnp+Onn35CUlISvvrqK2RnZ6OoqAjx8fEAgFWrVmH8+PEoLS1tcpTMniMiImqKFu05On3SMz1H8Z1dqu8777yD+fPno7i4GL169cLChQtx6623ulUPd0RERGD//v1NOvfV5UHEyMhIXH311c2qmCuOHTuGkpISZGRkyGlarRYDBw5Efn4+AKCgoAA1NTUOeeLj45GcnCzn2bFjB3Q6nRwYAUBaWhp0Op1DnuTkZDkwAoDMzExUV1ejoKDginWsrq6G0Wh0uIiIiHyJ7eBZdy9XrF69GlOnTsUzzzyDffv24dZbb8WQIUNw4sQJL7Wyca70BbkcHM2ePRuzZs3CxYsXXX3UJSUlJQCA2NhYh/TY2Fj5XklJCTQaDTp27NhgnpiYGKfyY2JiHPJc/j4dO3aERqOR89Rn3rx50Ol08pWQkOBiK4mIiLysFXbIXrBgASZMmICHHnoIPXr0wMKFC5GQkIB3333XS430LJfnHP3973/HL7/8gtjYWHTt2tVpQvbevXs9VjnAeWNJIUSj85wuz1Nf/ubkudzMmTMxbdo0+Wuj0cgAiYiIAtblIyRarRZardYhzWQyoaCgAE899ZRDekZGhjxi4+tcDo6Uq8m8Sa/XA7D26sTFxcnppaWlci+PXq+HyWRCWVmZQ+9RaWkp+vXrJ+c5c+aMU/lnz551KGfXrl0O98vKylBTU+PUo6RU3/8UREREvsS6CaR7mzjanr+8A2DWrFmYPXu2Q9rvv/8Os9nc4MiPr3M5OJo1a5Y36uGkW7du0Ov12Lx5M3r37g3AGo1u374dr7zyCgAgNTUVwcHB2Lx5M0aOHAkAKC4uxoEDB/Dqq68CANLT02EwGLB792784Q9/AADs2rULBoNBDqDS09Px0ksvobi4WA7ENm3aBK1Wi9TU1BZpLxERkTcIYb3cLQMAioqKHCZkN9RB0JyRH29y5b2bdfCsp1RUVODo0aPy18eOHcP+/fsRGRmJLl26YOrUqZg7dy6uvfZaXHvttZg7dy7atWuHMWPGAAB0Oh0mTJiA6dOnIyoqCpGRkZgxYwZSUlJwxx13ALBuDpWVlYXc3Fy8//77AKxL+bOzs5GUlATA2tXXs2dP5OTkYP78+Th//jxmzJiB3NxcrjojIiKqExER0ejvxejoaKjVaqdeIuXIT2vw+ITsyMhI/P77700utEuXLvjtt98azbdnzx707t1b7hmaNm0aevfujeeffx4A8OSTT2Lq1KmYOHEi+vTpg1OnTmHTpk3yHkcA8MYbb2DEiBEYOXIk+vfvj3bt2mH9+vXyHkcAsHLlSqSkpCAjIwMZGRm44YYb8NFHH8n31Wo1NmzYgJCQEPTv3x8jR47EiBEj8NprrzW5zURERL7IIoRHrqbSaDRITU3F5s2bHdI3b94sj9i0hq+++gpXXXVVk/I2aZ8jlUqFZcuWQafTNanQ+++/H4WFhU3aSyCQcJ8jIiJqipbc5+hksfu/k4xGIzrHNb2+q1evRk5ODt577z2kp6fjgw8+wKJFi/Djjz8iMTHRrbrU5+TJk/jyyy9x4sQJmEwmh3sLFixwubwmD6uNGzfO5cKJiIio7Rk1ahTOnTuHF154AcXFxUhOTsbGjRu9Ehh9/fXXGD58OLp164bDhw8jOTkZx48fhxACN998c7PKbFJwZLFYmlU4ERERtS6LsF7uluGqiRMnYuLEie69cRPMnDkT06dPxwsvvIDw8HCsWbMGMTExGDt2LLKysppVpu8es0tERERuE0J45PJVhw4dkke3goKCUFVVhfbt2+OFF16QV7e7isERERER+a2wsDBUV1cDsB4h9ssvv8j3XFlMptSqS/mJiIjIu1prWK2lpKWl4fvvv0fPnj0xdOhQTJ8+HYWFhVi7di3S0tKaVSaDIyIiogDnw7GN2xYsWICKigoA1vNfKyoqsHr1alxzzTV44403mlUmgyMiIqIAFug9R8ptg9q1a4d33nnH7TJdnnN02223Yc6cOU7pZWVluO2229yuEBEREVFTde/eHefOnXNKLy8vb/Z+iy73HG3btg2FhYXYt28fVq5cibCwMAD2c8+IiIjId3hitZkvr1Y7fvw4zGazU3p1dTVOnTrVrDKbNay2ZcsWPPLII0hLS8P69evRtWvXZr05EREReZel7nK3DF/z5Zdfyq//9a9/OZziYTab8fXXXzc7PmlWcBQXF4ft27fjwQcfxC233IJPP/0UPXr0aFYFiIiIiFw1YsQIAIAkSU6neAQHB6Nr1654/fXXm1W2y8GRJEkAAK1Wi5UrV+LFF19EVlYW/va3vzWrAkREROQ9Qlgvd8vwNbbTO7p164YffvgB0dHRHivb5eDo8nHHZ599Fj169ODZa0RERD4o0FerHTt2TH596dIlhISEuF2my6vVjh07hk6dOjmk3Xvvvdi1axc+/PBDtytERERE1FQWiwX//d//jauuugrt27fHr7/+CgB47rnnsHjx4maV6XJwlJiYKA+tKfXq1Yu9R0RERD4m0M9We/HFF7F06VK8+uqr0Gg0cnpKSgr+53/+p1ll8mw1IiKiAGbx0OWrli9fjg8++ABjx46FWq2W02+44Qb89NNPzSqTwRERERH5rVOnTuGaa65xSrdYLKipqWlWmQyOiIiIApiAfcVas6/WbkQDevXqhX//+99O6Z9++il69+7drDJ5thoREVEAswgBi5tzhtx93ptmzZqFnJwcnDp1ChaLBWvXrsXhw4exfPly/POf/2xWmew5IiIiCmDCQ5evGjZsGFavXo2NGzdCkiQ8//zzOHToENavX48777yzWWWy54iIiIj8WmZmJjIzMz1WHoMjIiKiABbom0DamEwmlJaWyjtn23Tp0sXlshgcERERBTIPHB/iy+NqR44cwYMPPoj8/HyHdCEEJEmC2Wx2uUwGR0REROS3xo8fj6CgIPzzn/9EXFxcvRtVu4rBERERUQCzQMDiZtePu8970/79+1FQUIDrr7/eY2UyOCIiIgpgwgPDaj68kh89e/bE77//7tEyuZSfiIiI/IrRaJSvV155BU8++SS2bduGc+fOOdwzGo3NKp89R0RERAEsEFerdejQwWFukRACt99+u0MeTsgmIiKiegXisNrWrVu9Wj6DIyIiIvIrAwcOlF+fOHECCQkJTqvUhBAoKipqVvmcc0RERBTAbKvV3L18Vbdu3XD27Fmn9PPnz6Nbt27NKpM9R0RERAEsEIfVlGxziy5XUVGBkJCQZpXJ4IiIiCiAWYSAxc3oxt3nvWHatGkAAEmS8Nxzz6Fdu3byPbPZjF27duGmm25qVtkMjoiIiMjv7Nu3D4C156iwsBAajUa+p9FocOONN2LGjBnNKpvBERERUQAzW6yXu2X4GtuKtQceeABvvvkmIiIiPFa2z0/I7tq1KyRJcroee+wxANYzVS6/l5aW5lBGdXU1Jk+ejOjoaISFhWH48OE4efKkQ56ysjLk5ORAp9NBp9MhJycH5eXlLdVMIiIir7ANq7l7+aolS5Z4NDAC/CA4+uGHH1BcXCxfmzdvBgDcd999cp6srCyHPBs3bnQoY+rUqVi3bh1WrVqF7777DhUVFcjOznbYGGrMmDHYv38/8vLykJeXh/379yMnJ6dlGklEREQ+w+eH1Tp16uTw9csvv4yrr77aYY8DrVYLvV5f7/MGgwGLFy/GRx99hDvuuAMAsGLFCiQkJGDLli3IzMzEoUOHkJeXh507d6Jv374AgEWLFiE9PR2HDx9GUlKSl1pHRETkXRYhYA7ACdne5PM9R0omkwkrVqzAgw8+6LBsb9u2bYiJicF1112H3NxclJaWyvcKCgpQU1ODjIwMOS0+Ph7JycnIz88HAOzYsQM6nU4OjAAgLS0NOp1OzlOf6upqj5zhQkRE5C3W40PcHVZr7Va0LL8Kjj7//HOUl5dj/PjxctqQIUOwcuVKfPPNN3j99dfxww8/4LbbbkN1dTUAoKSkBBqNBh07dnQoKzY2FiUlJXKemJgYp/eLiYmR89Rn3rx58hwlnU6HhIQED7SSiIiIWpPPD6spLV68GEOGDEF8fLycNmrUKPl1cnIy+vTpg8TERGzYsAH33HPPFcu6fNOo+jaQutLGUjYzZ86U91kArKcEM0AiIiJfEqir1bzJb4Kj3377DVu2bMHatWsbzBcXF4fExEQcOXIEAKDX62EymVBWVubQe1RaWop+/frJec6cOeNU1tmzZxEbG3vF99JqtdBqtc1pDhERUYsI1E0gvclvhtWWLFmCmJgYDB06tMF8586dQ1FREeLi4gAAqampCA4Olle5AUBxcTEOHDggB0fp6ekwGAzYvXu3nGfXrl0wGAxyHiIiImob/KLnyGKxYMmSJRg3bhyCguxVrqiowOzZs3HvvfciLi4Ox48fx9NPP43o6GjcfffdAACdTocJEyZg+vTpiIqKQmRkJGbMmIGUlBR59VqPHj2QlZWF3NxcvP/++wCAhx9+GNnZ2VypRkREfs3sgdVq7j7vb/wiONqyZQtOnDiBBx980CFdrVajsLAQy5cvR3l5OeLi4jB48GCsXr0a4eHhcr433ngDQUFBGDlyJKqqqnD77bdj6dKlUKvVcp6VK1diypQp8qq24cOH4+23326ZBhIREXmJBXB7tVkbm3IESYg2Fg56kdFohE6nw5mSEo/v1klERIHDaDQiVq+HwWDw2u8L2++kNXuOIqx9eOMPNKCy4gLu7XONV+vrS/xmzhERERFRS/CLYTUiIiJqHuGB1WptbZCJwREREVEAMwvr5W4ZbQmH1YiIiIgU2HNEREQUwLgJpOsYHBEREQUws0XA7OZafnef9zccViMiIqJW89JLL6Ffv35o164dOnTo0NrVAcDgiIiIKKDZhtXcvbzFZDLhvvvuw1/+8hevvYerOKxGREQUwHx9tdqcOXMAAEuXLvXem7iIPUdERERECuw5IiIiCmCeXK1mNBod0rVaLbRarVtl+yL2HBEREQUwi0V45AKAhIQE6HQ6+Zo3b1697zl79mxIktTgtWfPnpb8GFzCniMiIqIAZvHAnCPbSv6ioiKHg2ev1Gs0adIkjB49usEyu3bt6l6lvIjBERERETVJRESEQ3B0JdHR0YiOjm6BGnkHgyMiIqIA5us7ZJ84cQLnz5/HiRMnYDabsX//fgDANddcg/bt23vtfRvC4IiIiCiAmYWA2c3gxt3nG/L8889j2bJl8te9e/cGAGzduhWDBg3y2vs2hBOyiYiIqNUsXboUQginq7UCI4A9R0RERAFNudrMnTLaEgZHREREAcwMD+yQ7ZGa+A8OqxEREREpsOeIiIgogPn6ajVfxOCIiIgogPn6ajVfxGE1IiIiIgX2HBEREQUwi0XAzNVqLmFwREREFMDMHgiO3H3e3zA4IiIiCmAMjlzHOUdERERECuw5IiLXCUvz02ykev42a2oaETWZ2eJ+z4+5gW/lQMTgiIiIKIBxWM11/JOMiIiISIE9R0RUP9uQmHJorO61pEyzuDicVp/6hs5U1jShvGd7XV8aEdWLPUeuY3BEREQUwLjPkev4JxcRERGRAnuOiIiIAphZeGBYrY2drcbgiIganl9krnVOs5ib9GyTl/I3NJdIpXZKE+og53ych0RUL845cp1P/wSZPXs2JElyuPR6vXxfCIHZs2cjPj4eoaGhGDRoEH788UeHMqqrqzF58mRER0cjLCwMw4cPx8mTJx3ylJWVIScnBzqdDjqdDjk5OSgvL2+JJhIREZGP8engCAB69eqF4uJi+SosLJTvvfrqq1iwYAHefvtt/PDDD9Dr9bjzzjtx4cIFOc/UqVOxbt06rFq1Ct999x0qKiqQnZ0Ns9n+l++YMWOwf/9+5OXlIS8vD/v370dOTk6LtpOIiMgbbD1H7l5tic8PqwUFBTn0FtkIIbBw4UI888wzuOeeewAAy5YtQ2xsLD7++GM88sgjMBgMWLx4MT766CPccccdAIAVK1YgISEBW7ZsQWZmJg4dOoS8vDzs3LkTffv2BQAsWrQI6enpOHz4MJKSklqusURERB5WaxFQuxnc1DI48i1HjhxBfHw8tFot+vbti7lz56J79+44duwYSkpKkJGRIefVarUYOHAg8vPz8cgjj6CgoAA1NTUOeeLj45GcnIz8/HxkZmZix44d0Ol0cmAEAGlpadDpdMjPz28wOKqurkZ1dbX8tdFo9HDribzIYa8i67wih/2L6uYaSeYa53yWWqd8yvIk2+TNJs45EpLklAbbvCKV4seU7bWy3Lp8Dvsh2fJx7hER5xw1g0//5Ojbty+WL1+Of/3rX1i0aBFKSkrQr18/nDt3DiUlJQCA2NhYh2diY2PleyUlJdBoNOjYsWODeWJiYpzeOyYmRs5zJfPmzZPnKel0OiQkJDS7rUREROQbfLrnaMiQIfLrlJQUpKen4+qrr8ayZcuQlpYGAJBsf3HWEUI4pV3u8jz15W9KOTNnzsS0adPkr41GIwMkIiLyKdwE0nU+3XN0ubCwMKSkpODIkSPyPKTLe3dKS0vl3iS9Xg+TyYSysrIG85w5c8bpvc6ePevUK3U5rVaLiIgIh4uIiMiXmIXwyNWW+FVwVF1djUOHDiEuLg7dunWDXq/H5s2b5fsmkwnbt29Hv379AACpqakIDg52yFNcXIwDBw7IedLT02EwGLB79245z65du2AwGOQ8RAFFWKyXpVa+pFqT9aqpVlxVdddF+VJVV0JVXQlJcalMdVf1Bft1yQDVJQOki+VXvGx5VJcM9udM9vLk8hWXvS5ViquuvrY21JrsbbO11dWz3oioTfPpYbUZM2Zg2LBh6NKlC0pLS/Hiiy/CaDRi3LhxkCQJU6dOxdy5c3Httdfi2muvxdy5c9GuXTuMGTMGAKDT6TBhwgRMnz4dUVFRiIyMxIwZM5CSkiKvXuvRoweysrKQm5uL999/HwDw8MMPIzs7myvViIjI73FCtut8Ojg6efIk7r//fvz+++/o1KkT0tLSsHPnTiQmJgIAnnzySVRVVWHixIkoKytD3759sWnTJoSHh8tlvPHGGwgKCsLIkSNRVVWF22+/HUuXLoVabd91d+XKlZgyZYq8qm348OF4++23W7axREREXsDgyHWSEG1sINGLjEYjdDodzpSUcP4R+RbF0nt5Gb7iWBDbcn2p9pL9mVpTXVq1Ip81DTUmOU2YrM+IWpPiWWt5QrHZ6uUkxR8oCAq2pgVp7Pc1IdYXwfY0oba+FkFaxbO2tBBFPmt5UBwzIlT1bA1A1EqMRiNi9XoYDAav/b6w/U76r8XfQtOuvVtlmS5WYMWEAV6try/hTwkiIqIAxp4j1zE4IiIiCmBmYYHZ4t6iBHMbW9TgV6vViIiIiLyNPUdEgazurz2H4z5sc4kUx4LY5hpJNYr5RTUXrS+qq+Q0y6VKa7HV9rlJwpZWqzhm5PI5RxbF3COVda5R/XOOgu3vHxJm/Vdrn0ukktNCFU1s59BW62vrM0J5pInt5BHlkSI8XoTaAG4C6ToGR0RERAHMbBFQcc6RSxgcERERBbBaCyC5GdzUtq0pR5xzRERERKTEniOiQKOce9PEPY0kk3VekTzPCIC4WGEtotKoSLtgTauqtKfZ5hyZ7OVZTHVzjhpYISOp7H+bqTR1c4409vlFtjlHqtAw+3u1s76HKsy+z4rUru49GllNI881Us4zsu15xLlHFMA4rOY6BkdEREQBjMGR6/jnEhEREZECe46IAoVtWEl5VEhDy/ZN9iX6KpN1aEw5hGa5UG79t6LcnlZ331xxQU6ruWgtz3zJfnxIbZWprkrWOimPEbEt4VcOqwWFWo8AUYfYjwoJbmd9D3V7+1mJKtvQnWLbAFXdNgGqMPtftvUOsNUzdCZsPwGVR4pwiI0CDHuOXMfgiIiIKIBxnyPX8U8kIiIiIgX2HBEREQUws0W4vc8Rh9WIyD/ZjgpRHplRNx/HYdl+3REhymX78lwiwzl7mm3OkdGeZjJa5yZVl1fIaTWV1rlLtZX29zDXWOc9Wer+FWZFndTWDmtVsP3Hj7rudVCYfSl/cJj1iBBtB/tcJk3dnCPlHCahPJrE9h6SVPeinqNClGl1n5XymBHOOaJAI4SAcDO4EaJtBUf8KUBERESkwJ4jIiKiAGaxCLcnVLe1CdkMjoiIiAKYEMLtYbG2NqzG4IjInynnytheK44Kgbluvk6tfd6Oba6R7XgQQDG/qO5fwD7X6NI5g5xme20y2ucrmS5YX9dUVstpNVW2OUfO84FsVMFq+XVwqPVHUXCYVk7ThLezNuGSvVzbHCb7zCQ7SWUvT1X3WlLMHxK2+8o9jWxp9cxD4twjChTC4oE5R22s54jf/UREREQK7DkiIiIKYJxz5DoGR0T+TLlsv244zeGokLphNanWPjSFauvSe4ejQuqOCFEu27cNoV0sLbM/WreEX7mU/1KZdXm9qcL+vrWXrHUxm8x11bT/YJVU1mX2ao19GCwoxPqjSNPeXs+QuuNIbNsCAIDF7HwwiG2IzXYsibXAYOt71P1rvV93NIlaOawWZKuUnCTqW/JP5MeExXEEvrlltCX87iciIqJWcfz4cUyYMAHdunVDaGgorr76asyaNQsmk6nxh72IPUdEREQBzJdXq/3000+wWCx4//33cc011+DAgQPIzc1FZWUlXnvtNa+8Z1MwOCIiIgpgvjznKCsrC1lZWfLX3bt3x+HDh/Huu+8yOCKiZqpvKb/FPkdHqrXOA7LNPQIAyyXrESDi4gV7Wt38I9vxIIB9zpFyftHFUms+2zwjAKgqs85hqqm0zzmyL+W31kk5V0glHx9iH9W3LeU3VdrnCJlNFod/lWxlWMuxPqvV2OdQSRrrTCSL1r7oX6W1HkeCWo2cJtR1n5Wwv2+bm1xB5AKj0ejwtVarhVarvULu5jEYDIiMjPRoma7inCMiIqIAZtvnyN0LABISEqDT6eRr3rx5Hq3rL7/8grfeeguPPvqoR8t1FYMjIiKiQOaJwKguOCoqKoLBYJCvmTNn1vuWs2fPhiRJDV579uxxeOb06dPIysrCfffdh4ceesjrH0tDOKxGRERETRIREYGIiIhG802aNAmjR49uME/Xrl3l16dPn8bgwYORnp6ODz74wN1quo3BEZE/qpsXIyn3ObKY6/5VHB9ie11jn3Mkqq3zhSxV9vlF5grr/CPl/CLbESH17Wl08Xf78SFVtn2OLtrnHFXUWutlqvtr06yYy6m2bnMETd1+RwDQvu5ZTZW97sLsPAFUrbF2dquD7T+61CHW+Q5BIfa5RFJImPXfduH2h22fQZBiTpbt87Eojjmp+0yFcu4R9zwiP2YRApKbq80sLj4fHR2N6OjoJuU9deoUBg8ejNTUVCxZsgQqVet/vzE4IiIiCmBCeOBsNS8t5T99+jQGDRqELl264LXXXsPZs2fle3q93ivv2RQMjoiIiAKYLx88u2nTJhw9ehRHjx5F586dHd/TSwFZU7R+3xURERG1SePHj5c3qbz8ak3sOSLyZ5Z69jkyO8+pESb7vkTCts/RJfuco5qL1vs1lVVymumCdV5RfXsaVSnSKiqsc3kMNfa6XKr7K7OqnrPQbEIVexVV1c0v0tXz16mkts9Nss05Cgqxz3kKCgtxaAMAqOtpo+0zkLTt7Gm2zypIUU/bZ6o4qo3In1ksgOT2JpAeqoyfYHBEREQUwHz5+BBf5dPDavPmzcMtt9yC8PBwxMTEYMSIETh8+LBDnvHjxzvtnZCWluaQp7q6GpMnT0Z0dDTCwsIwfPhwnDx50iFPWVkZcnJy5I2tcnJyUF5e7u0mEhERkY/x6Z6j7du347HHHsMtt9yC2tpaPPPMM8jIyMDBgwcRFhYm58vKysKSJUvkrzUajUM5U6dOxfr167Fq1SpERUVh+vTpyM7ORkFBAdRqa9/5mDFjcPLkSeTl5QEAHn74YeTk5GD9+vUt0FIiF9V3xIUtTXnPNqxWq1jKX3ekiHKozXzJer+20p5WU1kNADBVKI4FqTsiRLls3zacZqy1v6/zUn77X51qyTpMVqVYqt8+yPnvNE3dewRX2n9M2epiq5uyzrY2KNtma6v1tfW+w1YH9X1ml98DuJSf/JqwuH8qTls7VcengyNboGKzZMkSxMTEoKCgAAMGDJDTtVrtFZf8GQwGLF68GB999BHuuOMOAMCKFSuQkJCALVu2IDMzE4cOHUJeXh527tyJvn37AgAWLVqE9PR0HD58GElJSV5qIRERkXdZLMIDc444rOazDAbrQZiXH0i3bds2xMTE4LrrrkNubi5KS0vlewUFBaipqUFGRoacFh8fj+TkZOTn5wMAduzYAZ1OJwdGAJCWlgadTifnISIiorbBp3uOlIQQmDZtGv74xz8iOTlZTh8yZAjuu+8+JCYm4tixY3juuedw2223oaCgAFqtFiUlJdBoNOjYsaNDebGxsSgpKQEAlJSUICYmxuk9Y2Ji5Dz1qa6uRnW1vXv/8tOKiYiIWpsv73Pkq/wmOJo0aRL+85//4LvvvnNIHzVqlPw6OTkZffr0QWJiIjZs2IB77rnniuUJISBJ9iXCytdXynO5efPmYc6cOa40g8izhPNSfuUxAbYjRaCYe2N7bTHZ02qrrPNxzDX2+Tg1dUd51F5yTquodV62r0yzLeGvqucIEMCaFqp2/t5SHiliK095pIitLjWKNFudbW1Qtk1dT7slxVEhts9K1PM5EgUKBkeu84thtcmTJ+PLL7/E1q1bnXbQvFxcXBwSExNx5MgRANbtx00mE8rKyhzylZaWIjY2Vs5z5swZp7LOnj0r56nPzJkzHU4nLioqcrVpRERE5GN8OjgSQmDSpElYu3YtvvnmG3Tr1q3RZ86dO4eioiLExcUBAFJTUxEcHIzNmzfLeYqLi3HgwAH069cPAJCeng6DwYDdu3fLeXbt2gWDwSDnqY9Wq5VPKG7qScVEREQtySKER662xKeH1R577DF8/PHH+OKLLxAeHi7P/9HpdAgNDUVFRQVmz56Ne++9F3FxcTh+/DiefvppREdH4+6775bzTpgwAdOnT0dUVBQiIyMxY8YMpKSkyKvXevTogaysLOTm5uL9998HYF3Kn52dzZVqRETk1zis5jqfDo7effddAMCgQYMc0pcsWYLx48dDrVajsLAQy5cvR3l5OeLi4jB48GCsXr0a4eHhcv433ngDQUFBGDlyJKqqqnD77bdj6dKl8h5HALBy5UpMmTJFXtU2fPhwvP32295vJFFzNHGfI2E2O/zrkKY4D8D22qKYc2SpseYzm8yKNMf9iwD7/CLHNOGQpqxtfd3Vasn5uJGwujlJFsWxJLa62OqmrHN97am33Yq0Ju9zROTHhPBAcMSeI9/R2H+M0NBQ/Otf/2q0nJCQELz11lt46623rpgnMjISK1ascLmOREREFFh8OjgiIiIi9wiLcHsTRw6rEVFgUyxlt7EPOTkPJSl/KFrq7te3Qt9cT09vfQNT9aXV96ztPSxm5XCZcz5bnR2Gy+wVrufdiNoWHjzrOp9erUZERETU0thzREREFMC4Ws11DI6IiIgCmMUiAB486xIGR0RtjUrtlCTVbWshqZ1H2iXFkR6quvv1nPwBtcNRO9YfpLbSGl/K71yg7T1Uijop62Kvu6ruX+d21ddWIqLGMDgiIiIKYMJihnBzcYK7z/sbBkdEREQBjMGR67hajYiIiEiBPUdE/kiq5+8aW5rinn0ukdo5TaWcy2N9rQq2/0hQBVvzqTVqRZo1n0Yx9ye0bs5PlWLzo9D6JiVdRpnHVl6oYn6RLc32nsq62OqmrHN97am33Yo0Uc9npiik0TYQ+QNhsXig56htHafD4IiIiCiACbO5/k1SXSyjLWFwREREFMCE8MCcI8HgiIj8hXLop+61UCyLF7al7EHB9nx1r1Uae1pQqAYAoFYMqwWHWl8HhTintb9YI6fZhtPaBzkPQ6kl6z3l8SC2ZfvKoTnbsyH1pNneU1kXZZqtzrY2OLStnnYLlXJYre796vkciajtYnBEREQUwLhazXUMjoiIiAIYgyPXsf+YiIiISIE9R0T+qIlL+aGqW+YeZJ+PI9XNvZE0IXKaOsR6PyjMnhYcpgUAaNpXy2mmSuuzmqpaOU1Xz5lLtvlEVeYrL/9VLtu3zTXSKZbta9oF19XDPm9I096WppXTbHW2tUHZNkkx50j+DFSKH3tcyk9tAHuOXMfgiIiIKIBxnyPX8U8jIiIiIgX2HBEREQUwi8UMuNlzZOGwGhH5DVU9+/Oo7d/WwjbnSDG/SAoJc/gXAILbXbD+GxYqp2nC2wEAQi6Z5DSzydq1Lsz1zDNS7H1UUWvNF1Z3RIgyu+3UkPr2ObLNMwKA0I4hdf/a6xRSl2arm7LOwe0abqPtMxDKOUe2z8phnhY71CmwcM6R6/hTgIiIiEiBPUdEREQBjD1HrmNwREREFMjMZgiVm8END54lIp8nn6OmnCtTd2aYYk6NZHsdrNgDSGude6MKtc/HUbcPBwBoOyjmF12y7m9krrHvaWSbc+RQlbpJRMGV9ve17YNkqbHmtyj2O1LV7W+kUuxpZDsrTbmnkW2ukW2ekbV+7a3lR7RzSrO1Qdk2W1uthTvvc2Sbf+R43loDex8R+SEh3J+Q3dYOnuV3PxEREZECe46IiIgCmLBY3O85amObQDI4IvJnUj1L+ZXDRrbjM8z2YTVV3fJ20e6SPc1kfa0x2dMsdcNplnqOAFFrVE6vTRX2pfy1l6zPmk3WH8hCccSIVLeEX62xD2UFhVjrbDseBLAPp9mGzZSvQ6J0cpomwtoeVViE/T3ahTu0FQCE2voZCMWRIvJnVd/nSBQghAf2OWprE7L5U4CIiIhIgT1HREREAcw6rObesBiH1YiIiChgcFjNdQyOiPyZYn6MsB2FIRR/4Vnq5hoF2ZfjS1rrEnnlHB3UWucLCcVeJopF8DLbMnx1sP1HR1DIRQBATWW1nFYjL+W/8g9UVbB9zpF9Kb9WTrMdEaJctm+ba6Scc6SKiLL+276DPc3WNq396BERZC3bNvfI+jq47l/Fj0LOOSJq8xgcERERBTD2HLmOwREREVEAs1jMkBgcuYT9x0REREQK7Dki8mf17c+jnD9Tz5wjYbHO4ZHaKY70sNj2I3L+61A590hVN9dIHWKfGxQUZs1RW2nfI8l25IhtrySh2CtJko8PsdfTNofJVhYABIdZ5wsp9zmS9zSqm2cEAKrwDg7/WttmfUYE2+crIajus1DMOZI/K+5zRAFMmC2A5GbPUT37nQUy/hS4zDvvvINu3bohJCQEqamp+Pe//93aVSIiImo2IcwQFjcvnq3Wdq1evRpTp07FM888g3379uHWW2/FkCFDcOLEidauGhERUbO4HRjVXd4yfPhwdOnSBSEhIYiLi0NOTg5Onz7ttfdrCgZHCgsWLMCECRPw0EMPoUePHli4cCESEhLw7rvvtnbViIiIAtLgwYPxv//7vzh8+DDWrFmDX375BX/+859btU6cc1THZDKhoKAATz31lEN6RkYG8vPz632muroa1dX2vV2MRqNX60hEROQqYTG7P+fIiz1HTzzxhPw6MTERTz31FEaMGIGamhoEBwc38KT3sOeozu+//w6z2YzY2FiH9NjYWJSUlNT7zLx586DT6eQrISGhJapKRETUZL4+rKZ0/vx5rFy5Ev369Wu1wAhgz5ETSZIcvhZCOKXZzJw5E9OmTZO/NhgM6NKlCy5cuODVOhLVy2JdGSZZFLth15is/9ZetKeZLtX9WymniUrra3OlPZ/lYlXdv/ZVaKYqa0+p6ZJJTquptr6ura6R08y1davUGlitJgkhp6kt1tdBQfa/14LV1h20q6vsvbOaIGuaKqhKTlOprKvP1Gr7SjdJWH+0CcXCNKGpq0uQ/Ye8CLaWLVSKH4Uq/lgk77P9nhCK7wOvMdfA7XcxW7+/Lx8h0Wq10Gq19T3hkr/97W94++23cfHiRaSlpeGf//yn22W6RZAQQojq6mqhVqvF2rVrHdKnTJkiBgwY0KQyioqKBABevHjx4sWrSVdRUZE3fqUJIYSoqqoSer3eY3Vt3769U9qsWbPqfe9Zs2Y1Wt4PP/wg5z979qw4fPiw2LRpk+jfv7+46667hMVi8dpn0xhJiJYIW/1D3759kZqainfeeUdO69mzJ/70pz9h3rx5jT5vsVhw+vRphIeHX7G36XJGoxEJCQkoKipCRERE4w8ECLab7W4L2G62+0qEELhw4QLi4+OhUnlvhsulS5dgMpkaz9gEop6RlCv1HP3+++/4/fffGyyva9euCAlxPsXx5MmTSEhIQH5+PtLT092rdDOx/1hh2rRpyMnJQZ8+fZCeno4PPvgAJ06cwKOPPtqk51UqFTp37tys946IiGhTP0Rs2O62he1uW9juhul0Oq/XJSQkpN4AxNuio6MRHR3drGdtfTbKBU8tjcGRwqhRo3Du3Dm88MILKC4uRnJyMjZu3IjExMTWrhoREVHA2b17N3bv3o0//vGP6NixI3799Vc8//zzuPrqq1ut1whgcORk4sSJmDhxYmtXg4iIKOCFhoZi7dq1mDVrFiorKxEXF4esrCysWrXKIxO9m4vBUSvTarWYNWtWq/5P0BrYbra7LWC72W5qWEpKCr755pvWroYTTsgmIiIiUuAmkEREREQKDI6IiIiIFBgcERERESkwOCIiIiJSYHDkIS+99BL69euHdu3aoUOHDvXmkSTJ6Xrvvfcc8hQWFmLgwIEIDQ3FVVddhRdeeMHp7J3t27cjNTUVISEh6N69u1MZALBmzRr07NkTWq0WPXv2xLp16zzWVqWmtPvEiRMYNmwYwsLCEB0djSlTpjjt2Opv7b5c165dnf7bPvXUUw55WvJz8DXvvPMOunXrhpCQEKSmpuLf//53a1epyWbPnu3031av18v3hRCYPXs24uPjERoaikGDBuHHH390KKO6uhqTJ09GdHQ0wsLCMHz4cJw8edIhT1lZGXJycuSDrHNyclBeXt4STQQAfPvttxg2bBji4+MhSRI+//xzh/st2c6mfK94SmPtHj9+vNN//7S0NIc8/thuakTrnFoSeJ5//nmxYMECMW3aNKHT6erNA0AsWbJEFBcXy9fFixfl+waDQcTGxorRo0eLwsJCsWbNGhEeHi5ee+01Oc+vv/4q2rVrJx5//HFx8OBBsWjRIhEcHCw+++wzOU9+fr5Qq9Vi7ty54tChQ2Lu3LkiKChI7Ny5s8XbXVtbK5KTk8XgwYPF3r17xebNm0V8fLyYNGmSX7f7comJieKFF15w+G974cKFVvkcfM2qVatEcHCwWLRokTh48KB4/PHHRVhYmPjtt99au2pNMmvWLNGrVy+H/7alpaXy/ZdfflmEh4eLNWvWiMLCQjFq1CgRFxcnjEajnOfRRx8VV111ldi8ebPYu3evGDx4sLjxxhtFbW2tnCcrK0skJyeL/Px8kZ+fL5KTk0V2dnaLtXPjxo3imWeeEWvWrBEAxLp16xzut1Q7m/K90pLtHjdunMjKynL473/u3DmHPP7YbmoYgyMPW7JkSYPB0eXfeErvvPOO0Ol04tKlS3LavHnzRHx8vHwA35NPPimuv/56h+ceeeQRkZaWJn89cuRIkZWV5ZAnMzNTjB492sXWNN2V2r1x40ahUqnEqVOn5LRPPvlEaLVaYTAYhBD+3W6bxMRE8cYbb1zxfkt+Dr7mD3/4g3j00Ucd0q6//nrx1FNPtVKNXDNr1ixx44031nvPYrEIvV4vXn75ZTnt0qVLQqfTiffee08IIUR5ebkIDg4Wq1atkvOcOnVKqFQqkZeXJ4QQ4uDBgwKAQyC/Y8cOAUD89NNPXmhVwy7/WdWS7WzK94q3XCk4+tOf/nTFZwKh3eSMw2otbNKkSYiOjsYtt9yC9957DxaLRb63Y8cODBw40GEDsczMTJw+fRrHjx+X82RkZDiUmZmZiT179qCmpqbBPPn5+V5q1ZXt2LEDycnJiI+Pd6hLdXU1CgoK5DyB0O5XXnkFUVFRuOmmm/DSSy85dIe35OfgS0wmEwoKCpzqnJGR0Sr/PzbXkSNHEB8fj27dumH06NH49ddfAQDHjh1DSUmJQ/u0Wi0GDhwot6+goAA1NTUOeeLj45GcnCzn2bFjB3Q6Hfr27SvnSUtLg06n84nPqSXb2ZTvlZa2bds2xMTE4LrrrkNubi5KS0vle4Hc7raMwVEL+u///m98+umn2LJlC0aPHo3p06dj7ty58v2SkhLExsY6PGP7uqSkpME8tbW18gnIV8pjK6Ml1VeXjh07QqPRNNom272G8vhKux9//HGsWrUKW7duxaRJk7Bw4UKHY2ha8nPwJb///jvMZrPP/P/YHH379sXy5cvxr3/9C4sWLUJJSQn69euHc+fOyW1oqH0lJSXQaDTo2LFjg3liYmKc3jsmJsYnPqeWbGdTvlda0pAhQ7By5Up88803eP311/HDDz/gtttukw9FDdR2t3UMjhpQ30TMy689e/Y0ubxnn30W6enpuOmmmzB9+nS88MILmD9/vkMeSZIcvhZ1k3GV6c3Nc3nalXi63fW97+X18YV2X86Vz+GJJ57AwIEDccMNN+Chhx7Ce++9h8WLF+PcuXNXrFt99fPU5+BrPPnfpaUNGTIE9957L1JSUnDHHXdgw4YNAIBly5bJeZrTvsb+2ze1nJbUUu30pc9i1KhRGDp0KJKTkzFs2DB89dVX+Pnnn+X/D67E39vd1vFstQZMmjQJo0ePbjBP165dm11+WloajEYjzpw5g9jYWOj1eqe/EGzdt7a/KK6UJygoCFFRUQ3mufyvkivxZLv1ej127drlkFZWVoaamppG2wS0bLsv587nYFvNcvToUURFRbXo5+BLoqOjoVarPfrfpbWFhYUhJSUFR44cwYgRIwBY/+qPi4uT8yjbp9frYTKZUFZW5tC7UFpain79+sl5zpw54/ReZ8+e9YnPybY6ryXa2ZTvldYUFxeHxMREHDlyBEDbaXdbw56jBkRHR+P6669v8AoJCWl2+fv27UNISIi8BD49PR3ffvutw1yVTZs2IT4+Xv4lnJ6ejs2bNzuUs2nTJvTp0wfBwcEN5rF9ozbGk+1OT0/HgQMHUFxc7FAXrVaL1NRUn2r35dz5HPbt2wcA8i+SlvwcfIlGo0FqaqpTnTdv3tzs/y6trbq6GocOHUJcXBy6desGvV7v0D6TyYTt27fL7UtNTUVwcLBDnuLiYhw4cEDOk56eDoPBgN27d8t5du3aBYPB4BOfU0u2synfK63p3LlzKCoqkr+320q725wWnf4dwH777Texb98+MWfOHNG+fXuxb98+sW/fPnk595dffik++OADUVhYKI4ePSoWLVokIiIixJQpU+QyysvLRWxsrLj//vtFYWGhWLt2rYiIiKh3KfcTTzwhDh48KBYvXuy0lPv7778XarVavPzyy+LQoUPi5Zdf9tqS9sbabVueevvtt4u9e/eKLVu2iM6dOzssT/XHdivl5+eLBQsWiH379olff/1VrF69WsTHx4vhw4fLeVryc/A1tqX8ixcvFgcPHhRTp04VYWFh4vjx461dtSaZPn262LZtm/j111/Fzp07RXZ2tggPD5fr//LLLwudTifWrl0rCgsLxf3331/vEvfOnTuLLVu2iL1794rbbrut3qXeN9xwg9ixY4fYsWOHSElJadGl/BcuXJC/fwHI/0/btlxoqXY25Xulpdp94cIFMX36dJGfny+OHTsmtm7dKtLT08VVV13l9+2mhjE48pBx48YJAE7X1q1bhRBCfPXVV+Kmm24S7du3F+3atRPJycli4cKFoqamxqGc//znP+LWW28VWq1W6PV6MXv2bHkZt822bdtE7969hUajEV27dhXvvvuuU30+/fRTkZSUJIKDg8X1118v1qxZ0yrtFsIaQA0dOlSEhoaKyMhIMWnSJIfl6v7YbqWCggLRt29fodPpREhIiEhKShKzZs0SlZWVDvla8nPwNf/4xz9EYmKi0Gg04uabbxbbt29v7So1mW0/n+DgYBEfHy/uuece8eOPP8r3LRaLmDVrltDr9UKr1YoBAwaIwsJChzKqqqrEpEmTRGRkpAgNDRXZ2dnixIkTDnnOnTsnxo4dK8LDw0V4eLgYO3asKCsra4kmCiGE2Lp1a73fy+PGjRNCtGw7m/K90hLtvnjxosjIyBCdOnUSwcHBokuXLmLcuHFObfLHdlPDJCEu236XiIiIqA3jnCMiIiIiBQZHRERERAoMjoiIiIgUGBwRERERKTA4IiIiIlJgcERERESkwOCIiIiISIHBERF5zPHjx+VDeW+66Sa3y7OVZTtih4ioJTA4IiKP27JlC77++mu3yykuLsbChQvdrxARkQsYHBGRx0VFRSEqKsrtcvR6PXQ6nQdqRETUdAyOiKheZ8+ehV6vx9y5c+W0Xbt2QaPRYNOmTS6VNX78eIwYMQJz585FbGwsOnTogDlz5qC2thZ//etfERkZic6dO+PDDz/0dDOIiFwW1NoVICLf1KlTJ3z44YcYMWIEMjIycP311+O//uu/MHHiRGRkZLhc3jfffIPOnTvj22+/xffff48JEyZgx44dGDBgAHbt2oXVq1fj0UcfxZ133omEhAQvtIiIqGnYc0REV3TXXXchNzcXY8eOxaOPPoqQkBC8/PLLzSorMjISf//735GUlIQHH3wQSUlJuHjxIp5++mlce+21mDlzJjQaDb7//nsPt4KIyDUMjoioQa+99hpqa2vxv//7v1i5ciVCQkKaVU6vXr2gUtl/5MTGxiIlJUX+Wq1WIyoqCqWlpW7XmYjIHQyOiKhBv/76K06fPg2LxYLffvut2eUEBwc7fC1JUr1pFoul2e9BROQJnHNERFdkMpkwduxYjBo1Ctdffz0mTJiAwsJCxMbGtnbViIi8hj1HRHRFzzzzDAwGA/7+97/jySefRI8ePTBhwoTWrhYRkVcxOCKiem3btg0LFy7ERx99hIiICKhUKnz00Uf47rvv8O6777Z29YiIvIbDakRUr0GDBqGmpsYhrUuXLigvL3e5rKVLlzqlbdu2zSnt+PHjLpdNRORpDI6IyOP69euHm266Cfn5+W6V0759e9TW1jZ7hRwRUXMwOCIij+ncuTOOHDkCANBqtW6Xt3//fgDWZf5ERC1FEkKI1q4EERERka/ghGwiIiIiBQZHRERERAoMjoiIiIgUGBwRERERKTA4IiIiIlJgcERERESkwOCIiIiISIHBEREREZECgyMiIiIihf8P4/i6uX03MLgAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "results[0]['theta_p'][0].plot.imshow(vmin=-3, vmax=3, cmap='RdBu_r')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Running the model forward in time\n", + "\n", + "Now, we can actually start integrating in time. At the first timestep, we don't have the data to use leapfrog integration, so we need a special first timestep that does an forward-in-time approach (also known as Euler forward differencing):" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "# Integrate the model, saving after desired timestep counts\n", + "model.take_first_timestep()\n", + "results.append(model.current_state())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "TODO: fix/adjust approach later\n", "\n", - "...TODO..." + "Now, the model (as written above) has been found to have a numerical instability that causes the model to \"explode\" within 30 seconds of simulated time. To showcase how such a model can go wrong, we save results every timestep, as follows:" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 20, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "for _ in range(3000):\n", + " model.integrate(1)\n", + " results.append(model.current_state())\n", + "\n", + "# Merge results\n", + "ds = xr.concat(results, 't')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "TODO: add holoviews visualization of \"what's going wrong\"" + ] }, { "cell_type": "markdown", @@ -213,12 +858,13 @@ "source": [ "## Summary\n", "\n", - "...TODO...\n", - "\n", - "Add one final `---` marking the end of your body of content, and then conclude with a brief single paragraph summarizing at a high level the key pieces that were learned and how they tied to your objectives. Look to reiterate what the most important takeaways were.\n", + "The preceding cookbook provided a crash-course on the development of a 2-D Mesoscale Numerical Model, starting with the basic dynamical equations through the conversion to computer code and ultimately model integration. This pre-configured \"test-case\" is intended to be used for educational purposes only. Though note that all the assumptions/configurations used herein may not be applicable for other situations. Future notebooks are to be included, demonstrating how various configurations (i.e., equation simplifications, discritization schemes, boundary conditions, grid styles, and spatiotemporal resolutions) all influences the performance and accuracy of the resultant output. Additionally, we plan to include explicit walk-throughs regarding stability analyses, corrections, and filtering techniques. Check-in regularly for updates to this Cookbook.\n", "\n", "### What's next?\n", - "Let Jupyter book tie this to the next (sequential) piece of content that people could move on to down below and in the sidebar. However, if this page uniquely enables your reader to tackle other nonsequential concepts throughout this book, or even external content, link to it here!" + "\n", + "TODO: reference later portions\n", + "\n", + "" ] }, { @@ -226,6 +872,10 @@ "metadata": {}, "source": [ "## Resources and references\n", + "\n", + "TODO: add these\n", + "\n", + "" ] } ], diff --git a/notebooks/intro/driver.py b/notebooks/intro/driver.py new file mode 100644 index 0000000..785e48e --- /dev/null +++ b/notebooks/intro/driver.py @@ -0,0 +1,218 @@ +import numpy as np +import xarray as xr + +from util import * + +metadata_attrs = { + 'u': { + 'units': 'm/s' + }, + 'w': { + 'units': 'm/s' + }, + 'theta_p': { + 'units': 'K' + }, + 'pi': { + 'units': 'dimensionless' + }, + 'x': { + 'units': 'm' + }, + 'x_stag': { + 'units': 'm' + }, + 'z': { + 'units': 'm' + }, + 'z_stag': { + 'units': 'm' + }, + 't': { + 'units': 's' + } +} + +class ModelDriver: + + coords = {} + prognostic_arrays = {} + base_state_arrays = {} + diagnostic_arrays = {} + params = {} + + def __init__(self, nx, nz, dx, dz, dt, **kwargs): + # Set parameters + self.nx = nx + self.nz = nz + self.dx = dx + self.dz = dz + self.dt = dt + for k, v in kwargs.items(): + if k.endswith('_tendency'): + setattr(self, k, v) + else: + self.params[k] = v + self.dtype = dtype = getattr(self, 'dtype', np.float32) + self.t_count = 0 + + # Define arrays + self.coords['x'] = np.arange(self.nx) * self.dx - self.nx * self.dx / 2 + self.coords['x_stag'] = np.arange(self.nx + 1) * self.dx - (self.nx + 1) * self.dx / 2 + self.coords['z'] = np.arange(self.nz) * self.dz + self.coords['z_stag'] = np.arange(self.nz + 1) * self.dz - self.dz / 2 + self.prognostic_arrays['u'] = np.zeros((3, nz, nx + 1), dtype=dtype) + self.prognostic_arrays['w'] = np.zeros((3, nz + 1, nx), dtype=dtype) + self.prognostic_arrays['theta_p'] = np.zeros((3, nz, nx), dtype=dtype) + self.prognostic_arrays['pi'] = np.zeros((3, nz, nx), dtype=dtype) + self.active_prognostic_variables = ['u', 'w', 'theta_p', 'pi'] + self.base_state_arrays['theta_base'] = np.zeros(nz, dtype=dtype) + self.base_state_arrays['PI_base'] = np.zeros(nz, dtype=dtype) + self.base_state_arrays['rho_base'] = np.zeros(nz, dtype=dtype) + ## Todo do we need others?? + + def initialize_isentropic_base_state(self, theta, pressure_surface): + # Set uniform potential temperature + self.base_state_arrays['theta_base'] = np.full( + self.base_state_arrays['theta_base'].shape, theta, dtype=self.dtype + ) + # Calculate pi based on hydrostatic balance (from surface) + self.base_state_arrays['PI_base'] = nondimensional_pressure_hydrostatic( + self.base_state_arrays['theta_base'], + self.coords['z'], + pressure_surface + ) + # Calculate density from theta and pi + self.base_state_arrays['rho_base'] = density_from_ideal_gas_law( + self.base_state_arrays['theta_base'], + self.base_state_arrays['PI_base'] + ) + + def initialize_warm_bubble(self, amplitude, x_radius, z_radius, z_center): + if np.min(self.base_state_arrays['theta_base']) <= 0.: + raise ValueError("Base state theta must be initialized as positive definite") + + # Create thermal bubble (2D) + theta_p, pi = create_thermal_bubble( + amplitude, self.coords['x'], self.coords['z'], x_radius, z_radius, 0.0, z_center, + self.base_state_arrays['theta_base'] + ) + # Ensure boundary conditions, and add time stacking (future, current, past) + self.prognostic_arrays['theta_p'] = np.stack([apply_periodic_lateral_zerograd_vertical(theta_p)] * 3) + self.prognostic_arrays['pi'] = np.stack([apply_periodic_lateral_zerograd_vertical(pi)] * 3) + + def prep_new_timestep(self): + for var in self.active_prognostic_variables: + # Future-current to current-past + self.prognostic_arrays[var][0:2] = self.prognostic_arrays[var][1:3] + + def take_first_timestep(self): + # check for needed parameters and methods + if not 'c_s_sqr' in self.params: + raise ValueError("Must set squared speed of sound prior to first timestep") + if not ( + getattr(self, 'u_tendency') + and getattr(self, 'w_tendency') + and getattr(self, 'theta_p_tendency') + and getattr(self, 'pi_tendency') + ): + raise ValueError("Must set tendency equations prior to first timestep") + + # Increment + self.t_count = 1 + + # Integrate forward-in-time + self.prognostic_arrays['u'][2] = ( + self.prognostic_arrays['u'][1] + + self.dt * apply_periodic_lateral_zerograd_vertical(self.u_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['pi'][1], self.base_state_arrays['theta_base'], self.dx, self.dz + )) + ) + self.prognostic_arrays['w'][2] = ( + self.prognostic_arrays['w'][1] + + self.dt * apply_periodic_lateral_zerow_vertical(self.w_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['pi'][1], self.prognostic_arrays['theta_p'][1], + self.base_state_arrays['theta_base'], self.dx, self.dz + )) + ) + self.prognostic_arrays['theta_p'][2] = ( + self.prognostic_arrays['theta_p'][1] + + self.dt * apply_periodic_lateral_zerograd_vertical(self.theta_p_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['theta_p'][1], self.base_state_arrays['theta_base'], self.dx, self.dz + )) + ) + self.prognostic_arrays['pi'][2] = ( + self.prognostic_arrays['pi'][1] + + self.dt * apply_periodic_lateral_zerograd_vertical(self.pi_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['pi'][1], self.base_state_arrays['theta_base'], + self.base_state_arrays['rho_base'], self.params['c_s_sqr'], self.dx, self.dz + )) + ) + + self.prep_new_timestep() + + def take_single_timestep(self): + # Check if initialized + if self.t_count == 0: + raise RuntimeError("Must run initial timestep!") + self.t_count += 1 + + # Integrate leapfrog + self.prognostic_arrays['u'][2] = ( + self.prognostic_arrays['u'][0] + + 2 * self.dt * apply_periodic_lateral_zerograd_vertical(self.u_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['pi'][1], self.base_state_arrays['theta_base'], self.dx, self.dz + )) + ) + self.prognostic_arrays['w'][2] = ( + self.prognostic_arrays['w'][0] + + 2 * self.dt * apply_periodic_lateral_zerow_vertical(self.w_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['pi'][1], self.prognostic_arrays['theta_p'][1], + self.base_state_arrays['theta_base'], self.dx, self.dz + )) + ) + self.prognostic_arrays['theta_p'][2] = ( + self.prognostic_arrays['theta_p'][0] + + 2 * self.dt * apply_periodic_lateral_zerograd_vertical(self.theta_p_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['theta_p'][1], self.base_state_arrays['theta_base'], self.dx, self.dz + )) + ) + self.prognostic_arrays['pi'][2] = ( + self.prognostic_arrays['pi'][0] + + 2 * self.dt * apply_periodic_lateral_zerograd_vertical(self.pi_tendency( + self.prognostic_arrays['u'][1], self.prognostic_arrays['w'][1], + self.prognostic_arrays['pi'][1], self.base_state_arrays['theta_base'], + self.base_state_arrays['rho_base'], self.params['c_s_sqr'], self.dx, self.dz + )) + ) + + self.prep_new_timestep() + + def integrate(self, n_steps): + for _ in range(n_steps): + self.take_single_timestep() + + def current_state(self): + """Export the prognostic variables, with coordinates, at current time.""" + data_vars = {} + for var in self.active_prognostic_variables: + if var == 'u': + dims = ('t', 'z', 'x_stag') + elif var == 'w': + dims = ('t', 'z_stag', 'x') + else: + dims = ('t', 'z', 'x') + data_vars[var] = xr.Variable(dims, self.prognostic_arrays[var][1:2].copy(), metadata_attrs[var]) + data_vars['x'] = xr.Variable('x', self.coords['x'], metadata_attrs['x']) + data_vars['x_stag'] = xr.Variable('x_stag', self.coords['x_stag'], metadata_attrs['x_stag']) + data_vars['z'] = xr.Variable('z', self.coords['z'], metadata_attrs['z']) + data_vars['z_stag'] = xr.Variable('z_stag', self.coords['z_stag'], metadata_attrs['z_stag']) + data_vars['t'] = xr.Variable('t', [self.t_count * self.dt], metadata_attrs['t']) + return xr.Dataset(data_vars) \ No newline at end of file diff --git a/notebooks/intro/util.py b/notebooks/intro/util.py new file mode 100644 index 0000000..8a585f6 --- /dev/null +++ b/notebooks/intro/util.py @@ -0,0 +1,70 @@ +# Helper functions for dry model +import numpy as np +import numba + +from constants import * + +@numba.njit() +def nondimensional_pressure_hydrostatic(theta, z, pressure_surface): + pi = np.zeros_like(theta) + # Start at (w level) surface as given + pi_sfc = (pressure_surface / p0)**(R_d / c_p) + # Go down, half level, for u level, using above-surface theta + pi[0] = pi_sfc + gravity * (z[1] - z[0]) / (2 * c_p * theta[1]) + # Now, integrate upward over full u levels + for i in range(1, pi.shape[0]): + theta_at_level = 0.5 * (theta[i] + theta[i - 1]) + pi[i] = pi[i - 1] - gravity * (z[i] - z[i - 1]) / (c_p * theta_at_level) + return pi + +@numba.njit() +def density_from_ideal_gas_law(theta, pi): + return p0 * pi ** (c_v / R_d) / (R_d * theta) + +@numba.njit() +def create_thermal_bubble(amplitude, x, z, x_radius, z_radius, x_center, z_center, theta_base): + # Coordinates in 2d + xx = np.broadcast_to(x[None, :], (z.shape[0], x.shape[0])) + zz = np.broadcast_to(z[:, None], (z.shape[0], x.shape[0])) + rad = np.sqrt(((zz - z_center) / z_radius)**2 + ((xx - x_center) / x_radius)**2) + # Create thermal bubble + theta_p = np.zeros_like(xx) + for k in range(rad.shape[0]): + for i in range(rad.shape[1]): + if rad[k, i] <= 1.0: + theta_p[k, i] = 0.5 * amplitude * (np.cos(np.pi * rad[k, i]) + 1.0) + # Create balanced pi, integrating downward from assumed zero at topmost level + pi = np.zeros_like(theta_p) + for k in range(rad.shape[0] - 2, -1, -1): + for i in range(rad.shape[1]): + integrand_trapz = 0.5 * ( + theta_p[k + 1, i] / theta_base[k + 1]**2 + + theta_p[k, i] / theta_base[k]**2 + ) + pi[k, i] = pi[k + 1, i] - gravity * (z[k + 1] - z[k]) / c_p * integrand_trapz + # Return results + return theta_p, pi + +@numba.njit() +def apply_periodic_lateral_zerograd_vertical(a): + # Bottom and top (no gradient) + for i in range(0, a.shape[1]): + a[0, i] = a[1, i] + a[a.shape[0] - 1, i] = a[a.shape[0] - 2, i] + # Left and right (mirrored) + for k in range(1, a.shape[0] - 1): + a[k, 0] = a[k, a.shape[1] - 2] + a[k, a.shape[1] - 1] = a[k, 1] + return a + +@numba.njit() +def apply_periodic_lateral_zerow_vertical(a): + # Bottom and top (fixed zero) + for i in range(0, a.shape[1]): + a[0, i] = a[1, i] = 0 + a[a.shape[0] - 1, i] = a[a.shape[0] - 2, i] = 0 + # Left and right (mirrored) + for k in range(1, a.shape[0] - 1): + a[k, 0] = a[k, a.shape[1] - 2] + a[k, a.shape[1] - 1] = a[k, 1] + return a