diff --git a/Project.toml b/Project.toml
index e8fd82d725..813c4508fd 100755
--- a/Project.toml
+++ b/Project.toml
@@ -101,7 +101,7 @@ StaticArrays = "1"
Statistics = "1.9"
StructArrays = "0.4, 0.5, 0.6, 0.7"
TimesDates = "0.3"
-XESMF = "0.1.4"
+XESMF = "0.1.5"
julia = "1.9"
oneAPI = "2.0.1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 54feb79795..9f12acdeaa 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,6 +13,7 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
TimesDates = "bdfc003b-8df8-5c39-adcd-3a9087f5df4a"
+XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5"
[compat]
CairoMakie = "0.11, 0.12, 0.13, 0.14"
diff --git a/docs/make.jl b/docs/make.jl
index d08c291b52..351a9a3394 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -11,15 +11,16 @@ Distributed.addprocs(2)
set_theme!(Theme(fontsize=20))
CairoMakie.activate!(type = "svg")
- using Oceananigans
using NCDatasets
+ using XESMF
+
+ using Oceananigans
+ using Oceananigans.AbstractOperations
using Oceananigans.Operators
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
- using Oceananigans.TurbulenceClosures
using Oceananigans.TimeSteppers
- using Oceananigans.AbstractOperations
-
+ using Oceananigans.TurbulenceClosures
using Oceananigans.BoundaryConditions: Flux, Value, Gradient, Open
bib_filepath = joinpath(dirname(@__FILE__), "oceananigans.bib")
@@ -173,18 +174,20 @@ format = Documenter.HTML(collapselevel = 1,
DocMeta.setdocmeta!(Oceananigans, :DocTestSetup, :(using Oceananigans); recursive=true)
-OceananigansNCDatasetsExt = if isdefined(Base, :get_extension)
- Base.get_extension(Oceananigans, :OceananigansNCDatasetsExt)
-else
- Oceananigans.OceananigansNCDatasetsExt
+modules = Module[]
+OceananigansNCDatasetsExt = isdefined(Base, :get_extension) ? Base.get_extension(Oceananigans, :OceananigansNCDatasetsExt) : Oceananigans.OceananigansNCDatasetsExt
+OceananigansXESMFExt = isdefined(Base, :get_extension) ? Base.get_extension(Oceananigans, :OceananigansXESMFExt) : Oceananigans.OceananigansXESMFExt
+
+for m in [Oceananigans, XESMF, OceananigansNCDatasetsExt, OceananigansXESMFExt]
+ if !isnothing(m)
+ push!(modules, m)
+ end
end
-makedocs(sitename = "Oceananigans.jl",
+makedocs(; sitename = "Oceananigans.jl",
authors = "Climate Modeling Alliance and contributors",
- format = format,
- pages = pages,
+ format, pages, modules,
plugins = [bib],
- modules = [Oceananigans, OceananigansNCDatasetsExt],
warnonly = [:cross_references],
draft = false, # set to true to speed things up
doctest = true, # set to false to speed things up
diff --git a/docs/src/appendix/library.md b/docs/src/appendix/library.md
index 1e62a6cd32..86beea13d8 100644
--- a/docs/src/appendix/library.md
+++ b/docs/src/appendix/library.md
@@ -71,6 +71,10 @@ Private = false
Modules = [Oceananigans.Fields]
Private = false
```
+```@docs
+XESMF.Regridder
+OceananigansXESMFExt.regrid!
+```
## Forcings
diff --git a/ext/OceananigansXESMFExt.jl b/ext/OceananigansXESMFExt.jl
index 031d11ffb1..72139bf705 100755
--- a/ext/OceananigansXESMFExt.jl
+++ b/ext/OceananigansXESMFExt.jl
@@ -9,24 +9,28 @@ using Oceananigans.Grids: AbstractGrid, λnodes, φnodes, Center, Face, total_le
import Oceananigans.Fields: regrid!
import XESMF: Regridder
-# permutedims below is used because Python's xESMF expects
-# 2D arrays with (x, y) coordinates with y varying in dim=1 and x varying in dim=2
+node_array(ξ::AbstractMatrix, Nx, Ny) = view(ξ, 1:Nx, 1:Ny)
-node_array(ξ::AbstractMatrix, Nx, Ny) = permutedims(view(ξ, 1:Nx, 1:Ny), (2, 1))
-vertex_array(ξ::AbstractMatrix, Nx, Ny) = permutedims(view(ξ, 1:Nx+1, 1:Ny+1), (2, 1))
-
-x_node_array(x::AbstractVector, Nx, Ny) = permutedims(repeat(view(x, 1:Nx), 1, Ny), (2, 1))
+x_node_array(x::AbstractVector, Nx, Ny) = repeat(view(x, 1:Nx), 1, Ny)
x_node_array(x::AbstractMatrix, Nx, Ny) = node_array(x, Nx, Ny)
-y_node_array(y::AbstractVector, Nx, Ny) = repeat(view(y, 1:Ny), 1, Nx)
+y_node_array(y::AbstractVector, Nx, Ny) = repeat(transpose(view(y, 1:Ny)), Nx, 1)
y_node_array(y::AbstractMatrix, Nx, Ny) = node_array(y, Nx, Ny)
-x_vertex_array(x::AbstractVector, Nx, Ny) = permutedims(repeat(view(x, 1:Nx+1), 1, Ny+1), (2, 1))
+vertex_array(ξ::AbstractMatrix, Nx, Ny) = view(ξ, 1:Nx+1, 1:Ny+1)
+
+x_vertex_array(x::AbstractVector, Nx, Ny) = repeat(view(x, 1:Nx+1), 1, Ny+1)
x_vertex_array(x::AbstractMatrix, Nx, Ny) = vertex_array(x, Nx, Ny)
-y_vertex_array(y::AbstractVector, Nx, Ny) = repeat(view(y, 1:Ny+1), 1, Nx+1)
+y_vertex_array(y::AbstractVector, Nx, Ny) = repeat(transpose(view(y, 1:Ny+1)), Nx+1, 1)
y_vertex_array(y::AbstractMatrix, Nx, Ny) = vertex_array(y, Nx, Ny)
+"""
+ xesmf_coordinates(grid::AbstractGrid, ℓx, ℓy, ℓz)
+
+Extract the coordinates (latitude/longitude) and the coordinates' bounds from
+`grid` at locations `ℓx, ℓy, ℓz`.
+"""
function xesmf_coordinates(grid::AbstractGrid, ℓx, ℓy, ℓz)
Nx, Ny, Nz = size(grid)
@@ -44,23 +48,26 @@ function xesmf_coordinates(grid::AbstractGrid, ℓx, ℓy, ℓz)
λv = x_vertex_array(λv, Nx, Ny)
φv = y_vertex_array(φv, Nx, Ny)
- return Dict("lat" => φ, # φ is latitude
- "lon" => λ, # λ is longitude
- "lat_b" => φv,
- "lon_b" => λv)
-end
-
-function xesmf_coordinates(dst_field::AbstractField, src_field::AbstractField)
-
- ℓx, ℓy, ℓz = Oceananigans.Fields.instantiated_location(src_field)
+ # Python's xESMF expects 2D arrays with (x, y) coordinates
+ # in which y varies in dim=1 and x varies in dim=2
+ # therefore we transpose the coordinate matrices
+ coords_dictionary = Dict("lat" => permutedims(φ, (2, 1)), # φ is latitude
+ "lon" => permutedims(λ, (2, 1)), # λ is longitude
+ "lat_b" => permutedims(φv, (2, 1)),
+ "lon_b" => permutedims(λv, (2, 1)))
- dst_grid = dst_field.grid
- src_grid = src_field.grid
+ return coords_dictionary
+end
- dst_coordinates = xesmf_coordinates(dst_grid, ℓx, ℓy, ℓz)
- src_coordinates = xesmf_coordinates(src_grid, ℓx, ℓy, ℓz)
+"""
+ xesmf_coordinates(field::AbstractField)
- return dst_coordinates, src_coordinates
+Extract the coordinates (latitude/longitude) and the coordinates' bounds from
+the `field`'s grid.
+"""
+function xesmf_coordinates(field::AbstractField)
+ ℓx, ℓy, ℓz = Oceananigans.Fields.instantiated_location(field)
+ return xesmf_coordinates(field.grid, ℓx, ℓy, ℓz)
end
"""
@@ -89,20 +96,24 @@ For more information, see the Python xESMF documentation at:
Example
=======
-```@example
+To create a regridder for two fields that live on different grids.
+
+```@example regridding
using Oceananigans
using XESMF
z = (-1, 0)
-tg = TripolarGrid(; size=(360, 170, 1), z, southernmost_latitude = -80)
-llg = LatitudeLongitudeGrid(; size=(360, 180, 1), z,
+tg = TripolarGrid(; size=(180, 85, 1), z, southernmost_latitude = -80)
+llg = LatitudeLongitudeGrid(; size=(170, 80, 1), z,
longitude=(0, 360), latitude=(-82, 90))
src_field = CenterField(tg)
dst_field = CenterField(llg)
-regridder = Oceananigans.Fields.Regridder(dst_field, src_field, method="conservative")
+regridder = XESMF.Regridder(dst_field, src_field, method="conservative")
```
+
+We can use the above regridder to regrid via [`regrid!`](@ref).
"""
function Regridder(dst_field::AbstractField, src_field::AbstractField; method="conservative")
@@ -117,7 +128,8 @@ function Regridder(dst_field::AbstractField, src_field::AbstractField; method="c
dst_Nz = size(dst_field)[3]
@assert src_field.grid.z.cᵃᵃᶠ[1:src_Nz+1] == dst_field.grid.z.cᵃᵃᶠ[1:dst_Nz+1]
- dst_coordinates, src_coordinates = xesmf_coordinates(dst_field, src_field)
+ dst_coordinates = xesmf_coordinates(dst_field)
+ src_coordinates = xesmf_coordinates(src_field)
periodic = Oceananigans.Grids.topology(src_field.grid, 1) === Periodic ? true : false
regridder = XESMF.Regridder(src_coordinates, dst_coordinates; method, periodic)
@@ -155,7 +167,7 @@ llg = LatitudeLongitudeGrid(; size=(360, 180, 1), z,
src_field = CenterField(tg)
dst_field = CenterField(llg)
-λ₀, φ₀ = 150, 30. # degrees
+λ₀, φ₀ = 150, 30 # degrees
width = 12 # degrees
set!(src_field, (λ, φ, z) -> exp(-((λ - λ₀)^2 + (φ - φ₀)^2) / 2width^2))
diff --git a/test/test_xesmf.jl b/test/test_xesmf.jl
index ea0c54e93e..ea1118878a 100644
--- a/test/test_xesmf.jl
+++ b/test/test_xesmf.jl
@@ -4,6 +4,8 @@ using XESMF
using SparseArrays
using LinearAlgebra
+gaussian_bump(λ, φ; λ₀=0, φ₀=0, width=10) = exp(-((λ - λ₀)^2 + (φ - φ₀)^2) / 2width^2)
+
for arch in archs
@testset "XESMF extension [$(typeof(arch))]" begin
@info "Testing XESMF regridding [$(typeof(arch))]..."
@@ -34,9 +36,9 @@ for arch in archs
src_field = CenterField(src_grid)
dst_field = CenterField(dst_grid)
- λ₀, φ₀ = 150, 30. # degrees
width = 12 # degrees
- set!(src_field, (λ, φ, z) -> exp(-((λ - λ₀)^2 + (φ - φ₀)^2) / 2width^2))
+ set!(src_field,
+ (λ, φ, z) -> gaussian_bump(λ, φ; λ₀=150, φ₀=30, width) - 2gaussian_bump(λ, φ; λ₀=270, φ₀=-20, width))
regridder = XESMF.Regridder(dst_field, src_field)
diff --git a/validation/regridding/latitude_longitude_regridding.jl b/validation/regridding/latitude_longitude_regridding.jl
index b2e3e0ab14..2b6ee441ad 100644
--- a/validation/regridding/latitude_longitude_regridding.jl
+++ b/validation/regridding/latitude_longitude_regridding.jl
@@ -56,4 +56,3 @@ scatter!(ax5, Array(interior(c_xy, 180, 60, :)), z)
scatter!(ax6, Array(interior(c_xyz, 180, 60, :)), z′)
display(fig)
-
diff --git a/validation/regridding/xesmf_python_regridding.ipynb b/validation/regridding/xesmf_python_regridding.ipynb
new file mode 100644
index 0000000000..0e25ae1db2
--- /dev/null
+++ b/validation/regridding/xesmf_python_regridding.ipynb
@@ -0,0 +1,1309 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "ed6fcf50-9e24-4396-a83c-f4ce2ff02978",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "'0.8.10'"
+ ]
+ },
+ "execution_count": 1,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "import numpy as np\n",
+ "import xesmf as xe\n",
+ "import xarray as xr\n",
+ "import matplotlib.pyplot as plt\n",
+ "xe.__version__"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "e216e11f-9e8c-4ff9-b8df-88dd672197f7",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "'8.8.1'"
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "import esmpy\n",
+ "esmpy.__version__"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "cc39740e-dc94-4cba-bd75-74768f11d602",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# 2-degree lat-lon grid\n",
+ "ds_in = xe.util.grid_global(2, 2, lon1=360)\n",
+ "\n",
+ "# 1-degree lat-lon grid\n",
+ "ds_out = xe.util.grid_global(1, 1, lon1=360)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "75b0dc72-7e4d-4abe-a8a5-4c7065252ff7",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "
<xarray.DataArray 'gaussian_field' (y: 90, x: 180)> Size: 130kB\n",
+ "array([[ 1.46996071e-55, 1.14821221e-54, 8.72318053e-54, ...,\n",
+ " -1.68476030e-18, -5.10267008e-19, -1.50311816e-19],\n",
+ " [ 7.56948141e-55, 5.91265531e-54, 4.49195360e-53, ...,\n",
+ " -4.33216641e-18, -1.31209264e-18, -3.86509465e-19],\n",
+ " [ 3.79107863e-54, 2.96127832e-53, 2.24973791e-52, ...,\n",
+ " -1.08344883e-17, -3.28146499e-18, -9.66636986e-19],\n",
+ " ...,\n",
+ " [ 9.12064795e-39, 7.12429883e-38, 5.41246159e-37, ...,\n",
+ " -6.03438672e-28, -1.82764780e-28, -5.38379037e-29],\n",
+ " [ 4.19026120e-39, 3.27308686e-38, 2.48662462e-37, ...,\n",
+ " -1.38437947e-28, -4.19290015e-29, -1.23512284e-29],\n",
+ " [ 1.87237470e-39, 1.46254487e-38, 1.11112239e-37, ...,\n",
+ " -3.08896814e-29, -9.35562483e-30, -2.75593158e-30]],\n",
+ " shape=(90, 180))\n",
+ "Coordinates:\n",
+ " lon (y, x) float64 130kB 1.0 3.0 5.0 7.0 ... 353.0 355.0 357.0 359.0\n",
+ " lat (y, x) float64 130kB -89.0 -89.0 -89.0 -89.0 ... 89.0 89.0 89.0\n",
+ "Dimensions without coordinates: y, x
1.47e-55 1.148e-54 8.723e-54 ... -3.089e-29 -9.356e-30 -2.756e-30
array([[ 1.46996071e-55, 1.14821221e-54, 8.72318053e-54, ...,\n",
+ " -1.68476030e-18, -5.10267008e-19, -1.50311816e-19],\n",
+ " [ 7.56948141e-55, 5.91265531e-54, 4.49195360e-53, ...,\n",
+ " -4.33216641e-18, -1.31209264e-18, -3.86509465e-19],\n",
+ " [ 3.79107863e-54, 2.96127832e-53, 2.24973791e-52, ...,\n",
+ " -1.08344883e-17, -3.28146499e-18, -9.66636986e-19],\n",
+ " ...,\n",
+ " [ 9.12064795e-39, 7.12429883e-38, 5.41246159e-37, ...,\n",
+ " -6.03438672e-28, -1.82764780e-28, -5.38379037e-29],\n",
+ " [ 4.19026120e-39, 3.27308686e-38, 2.48662462e-37, ...,\n",
+ " -1.38437947e-28, -4.19290015e-29, -1.23512284e-29],\n",
+ " [ 1.87237470e-39, 1.46254487e-38, 1.11112239e-37, ...,\n",
+ " -3.08896814e-29, -9.35562483e-30, -2.75593158e-30]],\n",
+ " shape=(90, 180))
"
+ ],
+ "text/plain": [
+ " Size: 130kB\n",
+ "array([[ 1.46996071e-55, 1.14821221e-54, 8.72318053e-54, ...,\n",
+ " -1.68476030e-18, -5.10267008e-19, -1.50311816e-19],\n",
+ " [ 7.56948141e-55, 5.91265531e-54, 4.49195360e-53, ...,\n",
+ " -4.33216641e-18, -1.31209264e-18, -3.86509465e-19],\n",
+ " [ 3.79107863e-54, 2.96127832e-53, 2.24973791e-52, ...,\n",
+ " -1.08344883e-17, -3.28146499e-18, -9.66636986e-19],\n",
+ " ...,\n",
+ " [ 9.12064795e-39, 7.12429883e-38, 5.41246159e-37, ...,\n",
+ " -6.03438672e-28, -1.82764780e-28, -5.38379037e-29],\n",
+ " [ 4.19026120e-39, 3.27308686e-38, 2.48662462e-37, ...,\n",
+ " -1.38437947e-28, -4.19290015e-29, -1.23512284e-29],\n",
+ " [ 1.87237470e-39, 1.46254487e-38, 1.11112239e-37, ...,\n",
+ " -3.08896814e-29, -9.35562483e-30, -2.75593158e-30]],\n",
+ " shape=(90, 180))\n",
+ "Coordinates:\n",
+ " lon (y, x) float64 130kB 1.0 3.0 5.0 7.0 ... 353.0 355.0 357.0 359.0\n",
+ " lat (y, x) float64 130kB -89.0 -89.0 -89.0 -89.0 ... 89.0 89.0 89.0\n",
+ "Dimensions without coordinates: y, x"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# 2D lon/lat\n",
+ "lon2d = ds_in['lon']\n",
+ "lat2d = ds_in['lat']\n",
+ "\n",
+ "def gaussian(ampl, λ0, φ0, width):\n",
+ " # Gaussian field\n",
+ " return xr.DataArray(\n",
+ " ampl * np.exp(-((lon2d - λ0)**2 + (lat2d - φ0)**2) / (2 * width**2)),\n",
+ " dims=('y', 'x'),\n",
+ " coords={'lon': lon2d, 'lat': lat2d},\n",
+ " name=\"gaussian_field\"\n",
+ " )\n",
+ "\n",
+ "# pick a center and width\n",
+ "λ0, φ0 = 150, 30 # center\n",
+ "width = 12 # degrees\n",
+ "\n",
+ "data_in = gaussian(1, 150, 30, width) + gaussian(-2, 270, -20, width)\n",
+ "data_in"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "1aab428f-baf2-49a5-abb6-89ddd87c3c41",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "regridder = xe.Regridder(ds_in, ds_out, \"conservative\", periodic=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "a707eb6a-3f6b-488c-a5d8-2c85abf33cc4",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "
<xarray.DataArray (y: 180, x: 360)> Size: 518kB\n",
+ "array([[ 1.46996071e-55, 1.46996071e-55, 1.14821221e-54, ...,\n",
+ " -5.10267008e-19, -1.50311816e-19, -1.50311816e-19],\n",
+ " [ 1.47119826e-55, 1.47119826e-55, 1.14917888e-54, ...,\n",
+ " -5.10429692e-19, -1.50359739e-19, -1.50359739e-19],\n",
+ " [ 7.56948141e-55, 7.56948141e-55, 5.91265531e-54, ...,\n",
+ " -1.31209264e-18, -3.86509465e-19, -3.86509465e-19],\n",
+ " ...,\n",
+ " [ 4.19026120e-39, 4.19026120e-39, 3.27308686e-38, ...,\n",
+ " -4.19290015e-29, -1.23512284e-29, -1.23512284e-29],\n",
+ " [ 1.87284498e-39, 1.87284498e-39, 1.46291221e-38, ...,\n",
+ " -9.36223373e-30, -2.75787839e-30, -2.75787839e-30],\n",
+ " [ 1.87237470e-39, 1.87237470e-39, 1.46254487e-38, ...,\n",
+ " -9.35562483e-30, -2.75593158e-30, -2.75593158e-30]],\n",
+ " shape=(180, 360))\n",
+ "Coordinates:\n",
+ " lon (y, x) float64 518kB 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5\n",
+ " lat (y, x) float64 518kB -89.5 -89.5 -89.5 -89.5 ... 89.5 89.5 89.5\n",
+ "Dimensions without coordinates: y, x\n",
+ "Attributes:\n",
+ " regrid_method: conservative
1.47e-55 1.47e-55 1.148e-54 ... -9.356e-30 -2.756e-30 -2.756e-30
array([[ 1.46996071e-55, 1.46996071e-55, 1.14821221e-54, ...,\n",
+ " -5.10267008e-19, -1.50311816e-19, -1.50311816e-19],\n",
+ " [ 1.47119826e-55, 1.47119826e-55, 1.14917888e-54, ...,\n",
+ " -5.10429692e-19, -1.50359739e-19, -1.50359739e-19],\n",
+ " [ 7.56948141e-55, 7.56948141e-55, 5.91265531e-54, ...,\n",
+ " -1.31209264e-18, -3.86509465e-19, -3.86509465e-19],\n",
+ " ...,\n",
+ " [ 4.19026120e-39, 4.19026120e-39, 3.27308686e-38, ...,\n",
+ " -4.19290015e-29, -1.23512284e-29, -1.23512284e-29],\n",
+ " [ 1.87284498e-39, 1.87284498e-39, 1.46291221e-38, ...,\n",
+ " -9.36223373e-30, -2.75787839e-30, -2.75787839e-30],\n",
+ " [ 1.87237470e-39, 1.87237470e-39, 1.46254487e-38, ...,\n",
+ " -9.35562483e-30, -2.75593158e-30, -2.75593158e-30]],\n",
+ " shape=(180, 360))
- regrid_method :
- conservative
"
+ ],
+ "text/plain": [
+ " Size: 518kB\n",
+ "array([[ 1.46996071e-55, 1.46996071e-55, 1.14821221e-54, ...,\n",
+ " -5.10267008e-19, -1.50311816e-19, -1.50311816e-19],\n",
+ " [ 1.47119826e-55, 1.47119826e-55, 1.14917888e-54, ...,\n",
+ " -5.10429692e-19, -1.50359739e-19, -1.50359739e-19],\n",
+ " [ 7.56948141e-55, 7.56948141e-55, 5.91265531e-54, ...,\n",
+ " -1.31209264e-18, -3.86509465e-19, -3.86509465e-19],\n",
+ " ...,\n",
+ " [ 4.19026120e-39, 4.19026120e-39, 3.27308686e-38, ...,\n",
+ " -4.19290015e-29, -1.23512284e-29, -1.23512284e-29],\n",
+ " [ 1.87284498e-39, 1.87284498e-39, 1.46291221e-38, ...,\n",
+ " -9.36223373e-30, -2.75787839e-30, -2.75787839e-30],\n",
+ " [ 1.87237470e-39, 1.87237470e-39, 1.46254487e-38, ...,\n",
+ " -9.35562483e-30, -2.75593158e-30, -2.75593158e-30]],\n",
+ " shape=(180, 360))\n",
+ "Coordinates:\n",
+ " lon (y, x) float64 518kB 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5\n",
+ " lat (y, x) float64 518kB -89.5 -89.5 -89.5 -89.5 ... 89.5 89.5 89.5\n",
+ "Dimensions without coordinates: y, x\n",
+ "Attributes:\n",
+ " regrid_method: conservative"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "data_out = regridder(data_in)\n",
+ "data_out"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "34c39ae3-33d2-4d39-99ee-180756ddcace",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9UAAAGGCAYAAAB1z6JvAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAfdZJREFUeJzt3Xl8E2XiP/DPJG3TFtpCOZpWClQFEYrAgiJ4ACJgOQTxxkUQ1h8uoiKeyKrFr1JBxbqgsPhVQBHh6wGr6wUooIiu5fJARVwrBWztgqUHlKZNnt8fyUwzSZom06STTD/v12teaebKkxTy6TPzHJIQQoCIiIiIiIiIgmbSuwBERERERERE0YqVaiIiIiIiIiKNWKkmIiIiIiIi0oiVaiIiIiIiIiKNWKkmIiIiIiIi0oiVaiIiIiIiIiKNWKkmIiIiIiIi0oiVaiIiIiIiIiKNWKkmIiIiIiIi0oiVaqImyM3NhSRJIT3n0KFDMXTo0JCek4iIqLmFIyPD4bfffkNubi727dund1GC8uuvv0KSJKxatarRfYP5XUiShNzc3KYVzs22bdsgSRK2bdsWsnMSRZoYvQtARGovvPCC3kUgIiJqMX777TfMnz8fXbt2Rd++ffUuTsDS09PxxRdf4KyzztK7KEQtHivVRBGmZ8+eeheBiIiIglBdXY34+PhmuTNvt9tRV1cHi8WCCy+8MOyvR0SNY/NvogC999576Nu3LywWC7KysvD000977fPGG29g4MCBSElJQWJiIs4880xMmzYtqNfxbP4tN+96+umnsXjxYmRlZaF169YYNGgQvvzyy6a+LSIioiYLJCOB0OSkw+HAokWL0KNHD1gsFnTs2BE333wzjhw5otqva9eumDp1qtfx7jm7bds2nH/++QCAW265BZIkNdr8edWqVZAkCZs2bcK0adPQoUMHJCYmoqamBgCwfv16DBo0CK1atULr1q0xatQo7N271+s8L774Irp37w6LxYKePXti7dq1mDp1Krp27arsI/8NsGjRIjz++OPIysqCxWLB1q1bG2z+HejvoqKiArfeeivatWuH1q1b44orrsBPP/3kc9+DBw9i0qRJ6NixIywWC84991w8//zzXvv9+OOPuOKKK5CYmIj27dvjtttuQ2VlZYOfJZFR8E41UQA+/vhjjB8/HoMGDcK6detgt9uxaNEi/P7778o+X3zxBa6//npcf/31yM3NRXx8PA4dOoRPPvkkJGV4/vnn0aNHD+Tn5wMAHn74YYwePRqFhYVISUkJyWsQEREFK5CMBEKXk3/961+xYsUKzJo1C2PHjsWvv/6Khx9+GNu2bcOePXvQvn37gM/1pz/9CStXrsQtt9yCv/3tbxgzZgwAoFOnTo0eO23aNIwZMwavvvoqTp48idjYWCxYsAB/+9vflPPZbDY89dRTuOSSS/DVV18prdFWrFiBGTNm4Oqrr8azzz6L8vJyzJ8/X6mYe/r73/+O7t274+mnn0ZycjK6devmc79AfxdCCEyYMAE7d+7EI488gvPPPx+ff/45cnJyvM75/fffY/DgwejcuTOeeeYZWK1WfPTRR7jzzjtx7NgxPProowCA33//HUOGDEFsbCxeeOEFpKWl4bXXXsOsWbMa/SyJop4gokYNHDhQZGRkiOrqamVdRUWFSE1NFfJ/o6effloAECdOnGjSaw0ZMkQMGTJEeV5YWCgAiN69e4u6ujpl/VdffSUAiNdff71Jr0dERNQUgWSkEKHJyR9++EEAEDNnzlSt//e//y0AiIceekhZ16VLFzFlyhSvc3jmbEFBgQAgVq5cGVAZVq5cKQCIm2++WbW+qKhIxMTEiDvuuEO1vrKyUlitVnHdddcJIYSw2+3CarWKgQMHqvY7dOiQiI2NFV26dFHWyX8DnHXWWcJms6n2l7e5lzvQ38UHH3wgAIjnnntOdc4nnnhCABCPPvqosm7UqFGiU6dOory8XLXvrFmzRHx8vPjjjz+EEEI88MADQpIksW/fPtV+I0aMEADE1q1bBZFRsfk3USNOnjyJgoICTJw4EfHx8cr6pKQkjBs3TnkuNx+77rrr8H//9384evRoSMsxZswYmM1m5fl5550HADh06FBIX4eIiChQgWYkEJqc3Lp1KwB4Neu+4IILcO655+Ljjz8O+pxaXX311arnH330Eerq6nDzzTejrq5OWeLj4zFkyBBl9OsDBw6gpKQE1113ner4zp0746KLLvL5WldeeSViY2P9lieY34X8Od50002q9ZMmTVI9P336ND7++GNcddVVSExMVL2v0aNH4/Tp00pXtK1bt6JXr17o06eP33MSGREr1USNKCsrg8PhgNVq9drmvu7SSy/Fxo0blUDt1KkTsrOz8frrr4ekHO3atVM9t1gsAJyDoxAREekh0IwEQpOTx48fB+Ac+dpTRkaGsr05eJZBbmJ9/vnnIzY2VrWsX78ex44dA1D/HtLS0rzO6Wudr9fyJZjfxfHjxxETE+P1t4Wv/erq6rBkyRKv9zR69GgAUL2vQF6byIjYp5qoEW3btoUkSSgpKfHa5rlu/PjxGD9+PGpqavDll18iLy8PkyZNQteuXTFo0KDmKjIREVGzCCYjgabnpFwJLC4u9ur3/Ntvv6n6U8fHx/vso3zs2LGg+l03xHOkb/mcb775Jrp06dLgcfJ78OznDPj+zHy9li/B/C7atWuHuro6HD9+XFWx9tyvbdu2MJvNmDx5Mm6//Xafr5uVlaWcM9B/B0RGwzvVRI1o1aoVLrjgArz99ts4ffq0sr6yshLvvvuuz2MsFguGDBmChQsXAoDPUT+JiIiinZaMBLTn5GWXXQYAWLNmjWp9QUEBfvjhBwwfPlxZ17VrV3zzzTeq/X766SccOHDAqyxA01t+jRo1CjExMfjPf/6DAQMG+FwA4JxzzoHVasX//d//qY4vKirCzp07Nb9+ML+LYcOGAQBee+011fq1a9eqnicmJmLYsGHYu3cvzjvvPJ/vSa6UDxs2DPv378fXX3/t95xERsQ71UQB+J//+R9cccUVGDFiBO655x7Y7XYsXLgQrVq1wh9//AEAeOSRR3DkyBEMHz4cnTp1wokTJ/Dcc88hNjYWQ4YM0fkdEBERhUcgGQmEJifPOecc/L//9/+wZMkSmEwm5OTkKKN/Z2Zm4u6771b2nTx5Mv785z9j5syZuPrqq3Ho0CEsWrQIHTp0UJ3zrLPOQkJCAl577TWce+65aN26NTIyMpCRkRHU59C1a1c89thjmDdvHn755RdcccUVaNu2LX7//Xd89dVXaNWqFebPnw+TyYT58+djxowZuOaaazBt2jScOHEC8+fPR3p6Okwm7fe8Av1djBw5Epdeeinuv/9+nDx5EgMGDMDnn3+OV1991euczz33HC6++GJccskl+Otf/4quXbuisrISP//8M959911l9PbZs2fj5ZdfxpgxY/D4448ro3//+OOPmt8PUdTQe6Q0omjxzjvviPPOO0/ExcWJzp07iyeffFI8+uijymia//rXv0ROTo4444wzRFxcnOjYsaMYPXq0+Oyzz4J6nYZG/37qqae89oXHCJ1ERER6aCwjhQhdTtrtdrFw4ULRvXt3ERsbK9q3by/+/Oc/i8OHD6v2czgcYtGiReLMM88U8fHxYsCAAeKTTz7xylkhhHj99ddFjx49RGxsbKPZKo/+XVBQ4HP7xo0bxbBhw0RycrKwWCyiS5cu4pprrhFbtmxR7bdixQpx9tlni7i4ONG9e3fx8ssvi/Hjx4t+/fop+/j7G8DX6N9CBPa7EEKIEydOiGnTpok2bdqIxMREMWLECPHjjz/6fP+FhYVi2rRp4owzzhCxsbGiQ4cOYvDgweLxxx9X7ff999+LESNGiPj4eJGamiqmT58u/vnPf3L0bzI8SQghdKrPExERERERgBMnTqB79+6YMGECVqxYoXdxiCgIbP5NRERERNSMSkpK8MQTT2DYsGFo164dDh06hGeffRaVlZW466679C4eEQWJlWqiZmK32+GvYYgkSap5qImIiFqSlpSTFosFv/76K2bOnIk//vgDiYmJuPDCC7F8+XL06tVL7+IRUZDY/JuomXTt2hWHDh1qcPuQIUOwbdu25isQERFRBGFOElG04p1qomby7rvv+pwvU5aUlNSMpSEiIooszEkiila8U01ERERERESkkfaJ8IiIiIiIiIhaOMM3/3Y4HPjtt9+QlJQESZL0Lg4RUcQQQqCyshIZGRkwmUJ3jfX06dOw2Wyajo2Li0N8fHzIykItA7OeiMi3cGR9U3IeMGbWG75S/dtvvyEzM1PvYhARRazDhw+jU6dOITnX6dOn0S6hNU7Brul4q9WKwsJCw4UthReznojIv1BlfVNzHjBm1hu+Ui0PavHzwYMc4IKIyE1lZSXO7tYtpN+NNpsNp2DHTTgDcUH2MLLBgddKjsJmsxkqaCn8mPVERL6FOuubkvOAcbPe8JVquRlYUlISkpOTdS4NEVHkCUdz2QSYECcFF7ZmDptJGjHriYj8C3XWa8l5wLhZb/hKNRERNT+zJMEcZICbIQEGDVsiIiIj0ZLzgHGznpVqIiIKOZMEmIPMWhNgyKAlIiIyGi05Dxg361mpJiKikNN8p5qIiIgiXpPuVBsQK9VERBRyZg1XsM3hKQoRERGFmJacB4yb9aGbmJSIiIiIiIioheGdaiIiCjk2/yYiIjIuNv9WY6WaiIhCjs2/iYiIjIvNv9VYqSYiopDjnWoiIiLj4p1qNVaqiUh3QsOXshaSMOAcDhFKQvCDdhgzZomICGDWG42WnJePMyJWqomIKOR4p5qIiMi4eKdajZVqIiIKOfapJiIiMi72qVZjpZqIwqq5mnsFIpCysNkYERFRcJj11NKxUk1ERCHnvIIdbPNvIiIiigZach4wbtazUk1ERCHH5t9ERETGxebfaqxUE1HIRVIzsGD5KzubiwWOA5URERlXNOc8wKwPBQ5UpsZKNRERhZxJwxVsLVNzEBERUfPTkvOAcbOelWoiIgo53qkmIiIyLt6pVjPqxQIiIiIiIiKisOOdaiIKmWjvY9WYht4f+19540BlRETGxKwngAOVeWKlmoiIQo6VaiIiIuNipVqNlWoiIgo59qkmIiIyLvapVmOlmoiazOhNwRrj/v7ZPMzJDA13qvnRERFFLGY9s96dlpwHjJv1rFQTEVHImTRcwTa18D/YiIiIooWWnJePMyKO/k1ERERERESkEe9UExFRyGkaqMyYF6+JiIgMR/NAZQbNel3vVNfV1eFvf/sbsrKykJCQgDPPPBOPPfYYHA6Hso8QArm5ucjIyEBCQgKGDh2K/fv361hqIqKGCUlSLS2VPIBJsAsREVGkY9Zrz3mjZr2uleqFCxdi+fLlWLp0KX744QcsWrQITz31FJYsWaLss2jRIixevBhLly5FQUEBrFYrRowYgcrKSh1LTkRE/shXsINdKLp8+umnGDduHDIyMiBJEjZu3Oi1zw8//IArr7wSKSkpSEpKwoUXXoiioiJle01NDe644w60b98erVq1wpVXXokjR44047sgIqJgac15o2a9rpXqL774AuPHj8eYMWPQtWtXXHPNNRg5ciR27doFwHmXOj8/H/PmzcPEiRORnZ2N1atX49SpU1i7dq2eRSciIj949bplOHnyJPr06YOlS5f63P6f//wHF198MXr06IFt27bh66+/xsMPP4z4+Hhln9mzZ2PDhg1Yt24dduzYgaqqKowdOxZ2u7253gYREQWJd6rVdO1TffHFF2P58uX46aef0L17d3z99dfYsWMH8vPzAQCFhYUoKSnByJEjlWMsFguGDBmCnTt3YsaMGTqVnIgATq8RCM/PqKVMw2GSpKBH+DTqiKBGlpOTg5ycnAa3z5s3D6NHj8aiRYuUdWeeeabyc3l5OV566SW8+uqruPzyywEAa9asQWZmJrZs2YJRo0aFr/BEFBBmfeNa4nRbWnJePs6IdL1T/cADD+DGG29Ejx49EBsbi379+mH27Nm48cYbAQAlJSUAgLS0NNVxaWlpyjZPNTU1qKioUC1ERETUvBwOB9577z10794do0aNQseOHTFw4EBVE/Hdu3ejtrZWdfE8IyMD2dnZ2Llzpw6lJiIiCp6uler169djzZo1WLt2Lfbs2YPVq1fj6aefxurVq1X7SR5XNIQQXutkeXl5SElJUZbMzMywlZ+IiHyTzJKmJRiB9Od1t23bNkiS5LX8+OOPTXin1JDS0lJUVVXhySefxBVXXIFNmzbhqquuwsSJE7F9+3YAzovncXFxaNu2repYfxfPAV5AJyLSm9acN2rW69r8+7777sODDz6IG264AQDQu3dvHDp0CHl5eZgyZQqsVisAZ+imp6crx5WWlnrdvZbNnTsXc+bMUZ5XVFSwYk1EEaOlNBEzmSWYTOFt/i33573llltw9dVXB3zcgQMHkJycrDzv0KFDUK9LgZFn8hg/fjzuvvtuAEDfvn2xc+dOLF++HEOGDGnwWH8XzwHnBfT58+eHtsBERCHQUrp9acl5wLhZr2ul+tSpUzCZ1DfLzWazEsRZWVmwWq3YvHkz+vXrBwCw2WzYvn07Fi5c6POcFosFFoslvAUnIiL/zCZIpiAbQ0nB/eHRWH/ehnTs2BFt2rQJ+jgKTvv27RETE4OePXuq1p977rnYsWMHAMBqtcJms6GsrEx1t7q0tBSDBw9u8Ny8gE5EpDMtOQ8YNut1bf49btw4PPHEE3jvvffw66+/YsOGDVi8eDGuuuoqAM5m37Nnz8aCBQuwYcMGfPfdd5g6dSoSExMxadIkPYtORER+SCYNTcI0XPHWol+/fkhPT8fw4cOxdevWZnnNliguLg7nn38+Dhw4oFr/008/oUuXLgCA/v37IzY2Fps3b1a2FxcX47vvvvNbqbZYLEhOTlYtRETUfDTlvIGzXtc71UuWLMHDDz+MmTNnorS0FBkZGZgxYwYeeeQRZZ/7778f1dXVmDlzJsrKyjBw4EBs2rQJSUlJOpaciIj8MZklmILsN2WCc3/P/rGhaoGUnp6OFStWoH///qipqcGrr76K4cOHY9u2bbj00kubfP6WqKqqCj///LPyvLCwEPv27UNqaio6d+6M++67D9dffz0uvfRSDBs2DB9++CHeffddbNu2DQCQkpKC6dOn45577kG7du2QmpqKe++9F71791ZGAyciosijJecB42a9JIRBG/q7VFRUICUlBb+XlPBKNlGIcZqN0NCrv1VFRQXSrFaUl5eH7PtR/s7d0K0PWpnNQR170m7HVQe/9lr/6KOPIjc31++xkiRhw4YNmDBhQlCvOW7cOEiShHfeeSeo48hp27ZtGDZsmNf6KVOmYNWqVQCAl19+GXl5eThy5AjOOecczJ8/H+PHj1f2PX36NO677z6sXbsW1dXVGD58OF544YWgmnMz64nCh1kfGkbJ+qbkPGDcrNf1TjURERmTZAq+r5X8B8fhw4dVwR/OcTIuvPBCrFmzJmznN7qhQ4eisWvz06ZNw7Rp0xrcHh8fjyVLlmDJkiWhLh4RNQEr0+SPlpwHjJv1rFQTEVHINaX5d3P2kd27d69qdgkiIiJqXFObfxst61mpJqKg8ep1aBlx+g0tc1FKCG7/xvrzzp07F0ePHsUrr7wCAMjPz0fXrl3Rq1cv2Gw2rFmzBm+99RbeeuutoF6XiIgoWEbLei05Dxg361mpJiKikHOGbZDNv+EIav9du3ap+vPKUyzJ/XmLi4tRVFSkbLfZbLj33ntx9OhRJCQkoFevXnjvvfcwevTooF6XiIiopdOS84Bxs54DlRFR0HinOrya6+p1OAcqe7//+WhlDu667Ul7HUbvLghpeahlYNYThR6zPryiNeubkvOAcbNe13mqiYiIiIiIiKIZm3+T8Yngmpk0SOI1KGoe8t2BaO5vJUkSJFOQfaodvCtCRBqEKucBZj01GyFJLS7nAeNmPSvVREQUciazCaYg+1qZBP+YJSIiigZach4wbtazUk1ERCGnafRvYcyr10REREajefRvg2Y9K9VkHAE0/5KCaCImPJuA+Tq2hTUT46AlzSuam4axUk1EYcGsDztmffOJ5u5erFSrsVJNREQhx+bfRERExsXm32rGfFdEREREREREzYB3qin6NdDMy2fzryCahDXUOEXVVEw+XwtrGkbNx7MZXtQ0EdPSLMygTcKIKASaMesbbBLOrKcwicqs19j826hZz0o1ERGFnEmSYApyqg0T+/ERERFFBS05Lx9nRKxUExFRyElmE6Qg+1pJDt4FIiIiigZach4wbtazUk1ERCFnMkswBdkszOQw5tVrIiIio9GS84Bxs56VaopOPvpLefWrcn/eUN8Uz2NU/aXtrnUe/VzcD5f393ceg+AUG5EhWqbf0DSllkGDlog0CLQPdbBZ79Vf2jvrPb+JWkrWM+cjRzRkveYptQya9cb6NiAiIiIiIiJqRrxTTUREIcc+1URERMbFPtVqrFRTdPFoeqVqBib/7NlURjj8Nxfz5NGcS2n2pTT/qm8q1mgTMYM1DaPIEelNw0xmaOhTHabCEFF0ccvoBvPb/bvPtS5kWe/RJNxv1jPnKYwiOeu15Dxg3KxnpZqIiEJOMkmQgpxqI9j9iYiISB9acl4+zohYqSYiopAzmUwwBdkszGTnHR8iIqJooCXnAeNmPSvVREQUcppG/9bQjIyIiIian+bRvw2a9bpeKujatSskSfJabr/9dgCAEAK5ubnIyMhAQkIChg4div379+tZZNKLcHj1sVL6TsnbhHAtzueSo05ZIC92m3ORn9fZ1IuPfSSPRf16fsoll40ojIQkcRoUIjIGt0z1maehyPoG9mk45xvJeo9tROHAnI98ulaqCwoKUFxcrCybN28GAFx77bUAgEWLFmHx4sVYunQpCgoKYLVaMWLECFRWVupZbCIiaoQ8KmiwC0WXTz/9FOPGjUNGRgYkScLGjRsb3HfGjBmQJAn5+fmq9TU1NbjjjjvQvn17tGrVCldeeSWOHDkS3oITEVGTaM15o2a9ru+qQ4cOsFqtyvKvf/0LZ511FoYMGQIhBPLz8zFv3jxMnDgR2dnZWL16NU6dOoW1a9fqWWwiImqEZDJpWii6nDx5En369MHSpUv97rdx40b8+9//RkZGhte22bNnY8OGDVi3bh127NiBqqoqjB07Fna7PVzFJiKiJtKa80bN+ojpU22z2bBmzRrMmTMHkiThl19+QUlJCUaOHKnsY7FYMGTIEOzcuRMzZszQsbTUbBqaQst9vcPue5ujzvXo8N7W0NQEdihTaMjTZChTZ7i+BCRTjLK9fpv6NHIjHSGZon56LTY5ig5CkiJqyg2TWcNAZQa9em1kOTk5yMnJ8bvP0aNHMWvWLHz00UcYM2aMalt5eTleeuklvPrqq7j88ssBAGvWrEFmZia2bNmCUaNGha3sFEF8Zb1nk2p/We9w+N7m6ztRvlbjJ+vdc169TX0qyX0bs57CLNKm19KS8/JxRhQxleqNGzfixIkTmDp1KgCgpKQEAJCWlqbaLy0tDYcOHWrwPDU1NaipqVGeV1RUhL6wRETkn5YmXgYN2pbM4XBg8uTJuO+++9CrVy+v7bt370Ztba3qAnpGRgays7Oxc+dOVqqJiCKV1qbcBs36iKlUv/TSS8jJyfFqGiZ5XDkTQnitc5eXl4f58+eHpYxERBQYyRR82Bq1SVhLtnDhQsTExODOO+/0ub2kpARxcXFo27atan1aWppycd0XXkAnItKXlpyXjzOiiKhUHzp0CFu2bMHbb7+trLNarQCcgZuenq6sLy0t9bp77W7u3LmYM2eO8ryiogKZmZlhKDU1J6mhkTUd9vptcnNvu/NRkp8L4fazZ9Mwj/NKpvqmW0pzL9d/E4fzYo4wuY4xxyhNxJTTmMxe5RZR2hSMok+kNQ2jlm337t147rnnsGfPHr8Xw33hBfSWyWfWezb79pX1ru+8gLLeM+N9ZL17zjt3aTjr5XIx66m5MOsjU0R8A6xcuRIdO3ZU9bXKysqC1WpVRgQHnP2ut2/fjsGDBzd4LovFguTkZNVCRETNi4OX0GeffYbS0lJ07twZMTExiImJwaFDh3DPPfega9euAJwX0G02G8rKylTHBnIBvby8XFkOHz4czrdCREQeOFCZmu53qh0OB1auXIkpU6YgJqa+OJIkYfbs2ViwYAG6deuGbt26YcGCBUhMTMSkSZN0LDERETXGOW2GufEdVcdwtGcjmTx5sjL4mGzUqFGYPHkybrnlFgBA//79ERsbi82bN+O6664DABQXF+O7777DokWLGjy3xWKBxWIJX+GJiMgvLTnvPM6YWa97pXrLli0oKirCtGnTvLbdf//9qK6uxsyZM1FWVoaBAwdi06ZNSEpK0qGkREQUKC1zURp17kojq6qqws8//6w8LywsxL59+5CamorOnTujXbt2qv1jY2NhtVpxzjnnAABSUlIwffp03HPPPWjXrh1SU1Nx7733onfv3l4VciIiihxa55w2atbr/q5GjhwJIQS6d+/utU2SJOTm5qK4uBinT5/G9u3bkZ2drUMpqdn5mk7Dfb3DXt+f2lHnXOzORbLXqpe60/WLrVq1mDwW1Xb34+pOe50X9jrltSXhmrbLVS6f5W/oPUUoIUmcYiNKRcLvzmQyaVqC8emnn2LcuHHIyMiAJEnYuHFjo8ds374d/fv3R3x8PM4880wsX75c4zskANi1axf69euHfv36AQDmzJmDfv364ZFHHgn4HM8++ywmTJiA6667DhdddBESExPx7rvvwqzhDghFGX9Z75bzjWa9Z177yXrPbf6yXn4tv1kf6HuKUJGQF6SN3r87rTlv1KzX/U41EREZT3PcqT558iT69OmDW265BVdffXWj+xcWFmL06NG49dZbsWbNGnz++eeYOXMmOnToENDx5G3o0KEQQQyW8+uvv3qti4+Px5IlS7BkyZIQloyIiMKpue5UR0vWs1JNRERRKScnBzk5OQHvv3z5cnTu3Bn5+fkAgHPPPRe7du3C008/zUo1ERFRBIqWrGelmiKPW7MpyWtaDKFeLzcFA5xNtQBIdpvqEXW1yjQbotY1r6lrig7h8BgswWRWpsqQYl2D4MjTbMTEOh/NcfVF9Si6Mu2GkK/COSA3zFGm25DLzuk3yMCacqfac87hUA1K9cUXX2DkyJGqdaNGjcJLL72E2tpaxMbGNvk1iCgAHs2jfWV9g1No+cr6Ote6QLJeng7LX9a75TzQWNY7y8msp5amqXeqjZb1/J9OREQhJ0kaptlw/fGZmZmJlJQUZcnLywtJmUpKSrymaUpLS0NdXR2OHTsWktcgIiJqCTTlvIGznneqiYgo5Jpyp/rw4cNITk5W1ody6iTJY1AXuT+w53oiIiJqWFPvVBst61mpJiKikGtKpTo5OVkVtKFitVpRUlKiWldaWoqYmBivqZ+IiIioYU2tVBst61mppuggjy7r2e/KNa0G4N2XWnL1qRK2ajhqqp37u/pdCfnRo0+1ZDJDcvWdFq5HyZLgfBQJDRev/gTq8sFU39GKqJnJU21IQYzOHComswmmIMM22P2DNWjQILz77ruqdZs2bcKAAQPYn5ooErhnvVvOAw2PmyLV1kDYnBkfSNZLcl/qkGW962dmPelEr6zXkvPyceGkV9azTzUREUWlqqoq7Nu3D/v27QPgnEZj3759KCoqAgDMnTsXN998s7L/bbfdhkOHDmHOnDn44Ycf8PLLL+Oll17Cvffeq0fxiYh0o/ccx0SBipas551qIiIKOckkQTIF2fzbFNwfeLt27cKwYcOU53PmzAEATJkyBatWrUJxcbESugCQlZWF999/H3fffTeef/55ZGRk4O9//zun0yIiIgqSlpyXjwtGtGQ9K9UUkbym1/Bc75o2Q3LUKVNo+Gr2DQCi+iQcp086f6457dp22nUe7ym1pLh453kszkeTw6MMqidysy/JVVyT66nzUUgm5T14TbcRoXjl2nj0aBrWlD7VgRo6dKgy+Igvq1at8lo3ZMgQ7NmzJ6jXIaLw8Jv1bjmvevSR9aLamfEBZb3c/NtP1nuloJ+s95xCK1qynoxHSFLE57x8XDCiJetZqSYiopBrjko1ERER6aO5KtXRgpVqIiIKOXn+ymCPISIiosinJefl44yIlWqKHB7Nv3xuU5p/y8+F0kRMGe3Trdk3ADhOVcBxstK1zdUkTB4h1KNpN0ym+hFAbfGufZzNxry+AiQTYHL9F5If5fOZ3Mormf2/J4N+uVDkac5m4JLZDJO5gX/7fo4hIoMLNOvdcx7wm/WOUxXOXQLJelclIKCsl/PZX9YrWc6sJ/1Fes7LxxkR/4cTERERERERacQ71UREFHLsU01ERGRc7FOtxko1ERGFHCvVRERExsVKtRor1RTZXH1CPKfdkJ+rptRyPTpcfajkqTUcJysh5Gk2XP2s7TXOqTgctXWqlzPFxsDs6q8l969SemK5puEwufpjmWIt9X285EfXPsKtnELIU3EE++abF6fSolCSTBoGKtMw4AkRGYCPrHfPeV+P7lmv9KUOIOtNsc4/ff1lvSrnAb9ZL5eTWU8tjZacl48zIlaqiYgo5HinmoiIyLh4p1qNlWoiIgo5ySQFX6k28Q4KERFRNNCS8/JxRsRKNUUcSTgannLDc2ot4QDszqZbotbZzEuZbqPGNaWG7bTSFKy2yvlYV13jOlz9OpLJhJgEZ/OuWHml3MzLdT7Exde/njnOuc7so1wNlF3+KhGcXoN0IiQp7NNtsPk3ETXEs0uXF/e/A+RHP1mvTKEVQNbL3zN+s9495wFmPUWd5phai82/1Yz5roiIiIiIiIiaAe9UExFRyEkmMyTXnZ9gjiEiIqLIpyXn5eOMiJVqIiIKPZNZaU4Z1DFEREQU+bTkvHycAene/Pvo0aP485//jHbt2iExMRF9+/bF7t27le1CCOTm5iIjIwMJCQkYOnQo9u/fr2OJSRee/ayFcC2O+sVhBxx2iLpa52Jz9rMSNdWw19TAXlODumrnYj9tg/20DbUnT6sW+2lb/T6uY0RNtXORz+c6Pxx29esLR3253MsdBTjFRssjJCm8v3eTSdtCRC2Xqi+1j5xvKOtdOR1I1svr/Ga9e843lvX+xoGJIGH/zqeIFNbfu9acN2jW6/quysrKcNFFFyE2NhYffPABvv/+ezzzzDNo06aNss+iRYuwePFiLF26FAUFBbBarRgxYgQqKyv1KzgREfklmc2aFooun376KcaNG4eMjAxIkoSNGzcq22pra/HAAw+gd+/eaNWqFTIyMnDzzTfjt99+U52jpqYGd9xxB9q3b49WrVrhyiuvxJEjR5r5nRARUTC05rxRs17XSvXChQuRmZmJlStX4oILLkDXrl0xfPhwnHXWWQCcd6nz8/Mxb948TJw4EdnZ2Vi9ejVOnTqFtWvX6ll0IiLyR24WFuxCUeXkyZPo06cPli5d6rXt1KlT2LNnDx5++GHs2bMHb7/9Nn766SdceeWVqv1mz56NDRs2YN26ddixYweqqqowduxY2F2jPRMRUQTSmvMGzXpd+1S/8847GDVqFK699lps374dZ5xxBmbOnIlbb70VAFBYWIiSkhKMHDlSOcZisWDIkCHYuXMnZsyYoVfRqbl4NqnymAILwqFMzSEcdtUjlEcHHLV1rm3Ofe2u5w5bnddLmuKc/y3kY8zya/o4v/LaDZXT5PYeJGN+iRBRy5WTk4OcnByf21JSUrB582bVuiVLluCCCy5AUVEROnfujPLycrz00kt49dVXcfnllwMA1qxZg8zMTGzZsgWjRo0K+3ugCOAv6+UpqvxlvWv/hrJe2OvPZ4dznb+s9/w7wm/Wy7enmPVELZqud6p/+eUXLFu2DN26dcNHH32E2267DXfeeSdeeeUVAEBJSQkAIC0tTXVcWlqass1TTU0NKioqVAsRETUzk0nD1Wtj9rOieuXl5ZAkSenmtXv3btTW1qounmdkZCA7Oxs7d+7UqZRERNQoTTlv3KzX9U61w+HAgAEDsGDBAgBAv379sH//fixbtgw333yzsp/k0cFeCOG1TpaXl4f58+eHr9BERNQoyWSCFGRwBrs/RZfTp0/jwQcfxKRJk5CcnAzAefE8Li4Obdu2Ve3r7+I54LyAXlNTozznBXQioualJefl44xI13eVnp6Onj17qtade+65KCoqAgBYrVYA8ArW0tJSr7vXsrlz56K8vFxZDh8+HIaSExGRX5KGq9dsNmlYtbW1uOGGG+BwOPDCCy80ur+/i+eA8wJ6SkqKsmRmZoayuERE1BgtOW/grNe1Un3RRRfhwIEDqnU//fQTunTpAgDIysqC1WpV9cmy2WzYvn07Bg8e7POcFosFycnJqoWIiJoZBy8hl9raWlx33XUoLCzE5s2bVblstVphs9lQVlamOsbfxXOAF9CJiHTHgcpUdK1U33333fjyyy+xYMEC/Pzzz1i7di1WrFiB22+/HYCz2ffs2bOxYMECbNiwAd999x2mTp2KxMRETJo0Sc+iExGRH3KzsGAXMha5Qn3w4EFs2bIF7dq1U23v378/YmNjVRfPi4uL8d133zV48RzgBXQiIr1pzXmjZr2ufarPP/98bNiwAXPnzsVjjz2GrKws5Ofn46abblL2uf/++1FdXY2ZM2eirKwMAwcOxKZNm5CUlKRjyYmIiKiqqgo///yz8rywsBD79u1DamoqMjIycM0112DPnj3417/+BbvdrnTnSk1NRVxcHFJSUjB9+nTcc889aNeuHVJTU3Hvvfeid+/eymjgREREkU7XSjUAjB07FmPHjm1wuyRJyM3NRW5ubvMViiKH5HE1S766Za/fLuR9XM1JJLlZifJogik2xrXNua851vc/fXNsjLKPfIzymj7Or7x2Q+X0tS2CCD99FomaREsTL4M2CTOyXbt2YdiwYcrzOXPmAACmTJmC3NxcvPPOOwCAvn37qo7bunUrhg4dCgB49tlnERMTg+uuuw7V1dUYPnw4Vq1aBbOZ/x5aDF8Zaldv85v1HrntmfXyNFru6/xlveffEcx6Ih+0NuU2aNbrXqkmIiIDkqfaCPYYA5ArloFYvHhxGEsSfkOHDoUQosHt/rbJ4uPjsWTJEixZsiSURSMionDSkvPycQbESjUREYWcZDZDCvJOY7D7R6q9e/eqnu/evRt2ux3nnHMOAOeAnGazGf3799ejeERERE2mJefl44yIlWqKDp7NquSmTJKpfpvcJCwm1vkYF+98tCTAXFcLAIhJqFOdxhSn/i8gmUyISbAAAMwWi3K86nyu8zunBfBoEubZxCqCm4MRhZXJFPzVaINcvd66davy8+LFi5GUlITVq1crczGXlZXhlltuwSWXXKJXEYkik3tmuue8+6OPrJdzOpCsl5t9+8t6Vc77KoN71jPnqaXSkvPycQZkzHdFRET6asZpNl544QVkZWUhPj4e/fv3x2effdbgvtu2bYMkSV7Ljz/+qPWd+vXMM88gLy9PqVADQNu2bfH444/jmWeeCctrEhERhV0zTqkVyTkvY6WaiIii1vr16zF79mzMmzcPe/fuxSWXXIKcnBwUFRX5Pe7AgQMoLi5Wlm7duoWlfBUVFfj999+91peWlqKysjIsr0lERGQUkZ7zMlaqiYgo5CTXCLrBLsFavHgxpk+fjr/85S8499xzkZ+fj8zMTCxbtszvcR07doTValWWcI00fdVVV+GWW27Bm2++iSNHjuDIkSN48803MX36dEycODEsr0lERBRuWnM+2KyP9JyXsU81RRwhmaD0VhJ29Uavfk0mwPWfRIp19osScj8ri6sPtC0ecDjP4+olpUyh4aj16HcVG1PfvyqhlfNR6ZvtfIR8/lgLhPwf1Fe5PHlOC6IjTq9BYSdp6GsV5P8Nm82G3bt348EHH1StHzlyJHbu3On32H79+uH06dPo2bMn/va3v6mmhQql5cuX495778Wf//xn1Na6+nvGxGD69Ol46qmnwvKaRJFOzkG/We+Zp36yXrK58jmArJfX+c16t5wHwKwn8kVLzsvHBSgacl7GSjUREYWclqvR8v4VFRWq9RaLBRbXH8Dujh07BrvdjrS0NNX6tLQ0lJSU+HyN9PR0rFixAv3790dNTQ1effVVDB8+HNu2bcOll14aVHkDkZiYiBdeeAFPPfUU/vOf/0AIgbPPPhutWrUK+WsRERE1F60tzILJ+mjIeRkr1UREFHpNmKc6MzNTtfrRRx9Fbm5ug4dJHndjhBBe62TnnHOOMrUVAAwaNAiHDx/G008/HdawbdWqFc4777ywnZ+IiKhZNXGe6mCyPhpynpVqimyu/zBKUzHPZlWmGOci/4z6aTFMDodzvcMOh3w+139+c0y189GhbHFtN3lNoWVqleR8jG+lOr8wxUC4XlN59GgSJiST9zRbRC1BE6bUOnz4MJKTk5XVvu5SA0D79u1hNpu9rlaXlpZ6XdX258ILL8SaNWuCK6sfwfSVfvvtt0P2ukRRy0fWq3Lex6Mq613NvgPKetf3jL+sd895n49uWV//M7OeWpgmTqkVSNZHas77wko1ERFFlOTkZFXQNiQuLg79+/fH5s2bcdVVVynrN2/ejPHjxwf8env37kV6erqmsvqSkpISsnMREREZUSBZH6k57wsr1UREFHKS2QwpyJE2g90fAObMmYPJkydjwIABGDRoEFasWIGioiLcdtttAIC5c+fi6NGjeOWVVwAA+fn56Nq1K3r16gWbzYY1a9bgrbfewltvvRX0azdk5cqVITsXERFRJNKS8/JxwYjEnPeFlWqKHHITKuFoeJvSHMz16JDqm4TJo3WKBNWhJkBpCiZqTjsf49QjhdbvbPYa7Vtp9q2MEOpqFhYTW//acrMvk0c5/Y1wGAEjgxKFjcmsoU918OF8/fXX4/jx43jsscdQXFyM7OxsvP/+++jSpQsAoLi4WDWXpc1mw7333oujR48iISEBvXr1wnvvvYfRo0cH/dqBqqurw7Zt2/Cf//wHkyZNQlJSEn777TckJyejdevWYXtdoogkmXznvLxNfnTPecBv1itpGkjWu/bxl/WqnHd/bV9Z31iWM+vJqLTkvHxcEKIh5wFAEkKIsL6CzioqKpCSkoLfS0oCak5IEUA4IMmBqzw6/5lKDte0GHab83mdDZK91vWzM0Sl2hrnITZnXypRfRKO0yedP8tBa3M+NqlSHWuBiHHuI8zO4BUxcc7zmJ2PwhRTH8Ke/cF1DFpOs0GA8/vRmpaG8vLykH0/yt+5fxR8gOTWwY1wXVF1Eqnn54S0PHo7dOgQrrjiChQVFaGmpgY//fQTzjzzTMyePRunT5/G8uXL9S6iITDro4wr2/1mvVvOA/Cb9aLamfEBZX0wlWpl+q6Gs154VriZ9RRhQp31Tcl5wJhZD7hd3CMiIgoVyWTStBjNXXfdhQEDBqCsrAwJCfV31q666ip8/PHHOpaMiIhIO605b8SsB9j8m4iIwkHS0CxM0tCMLMLt2LEDn3/+OeLi4lTru3TpgqNHj+pUKiIioibSkvPycQbESjVFJGVaDWWFXb1emd7CoTQXk5tcy9wbPZnkq2Ku5l6iztmMTHg0/5ZMZkhyHyq535Yy7YZHczC3Zl/CY6oP4d7cy7MpmI6itSlYMJ1UovQtkkE5HA7Y7Xav9UeOHEFSUpIOJSKKHP6yXpXzQMiyXpKbf/vJevecd5aBWR9uzHmKdvr/zyciIuORpPpBfAJejPeX0ogRI5Cfn688lyQJVVVVePTRR8M+aAoREVHYaMp5Y2Y9wDvVREQUDoGMiuvrGIN59tlnMWzYMPTs2ROnT5/GpEmTcPDgQbRv3x6vv/663sUjIiLSRkvOy8cZECvVFB2Uq1oeo2maY9BoiyHJBJPclMs1WqgkNwXzMfq3MjKofIzczMvVRExpDmaOU0YChVke/dNzSi1jXo0LtVDOQeDvXPx1NB8hmYJuBhkJzSZDLSMjA/v27cPrr7+OPXv2wOFwYPr06bjppptUA5cREdRZ75bzAPxnvWvfgLJe7gPqJ+vdc9756CfrGSwBC1XWN3Ye/kqah5acl48zIlaqiYgo9HinWpGQkIBp06Zh2rRpeheFiIgoNHinWoWVaiIiCj1JCv52gUFuL7zzzjvIyclBbGws3nnnHb/7Xnnllc1UKiIiohDSkvPycQbESjUREVEITZgwASUlJejYsSMmTJjQ4H6SJPkcGZyIiIiii66V6tzcXMyfP1+1Li0tDSUlJQAAIQTmz5+PFStWoKysDAMHDsTzzz+PXr166VFcai6SCRDOKTS8ptuAx3pT/T9hry42cvMSUwzgqHP+LE/F4Tq/5HpUzuHWlEWYnf2tvKbQkB/NsfX9qzyn15Cvwvnqb6JDs5dInF7DX58oR4g6Xpk83rfnaSPwYzEOk8m5BHuMATgcDp8/E5GLnIN+st5XzgMNZL28TwBZLzzGPvGV9aqcB/xnfUNTaTHrAeif9RH4kRiHlpyXjzMg3d9Vr169UFxcrCzffvutsm3RokVYvHgxli5dioKCAlitVowYMQKVlZU6lpiIiBojD2AS7GIEqampOHbsGABg2rRphs6sTz/9FOPGjUNGRgYkScLGjRtV24UQyM3NRUZGBhISEjB06FDs379ftU9NTQ3uuOMOtG/fHq1atcKVV16JI0eONOO7ICKiYGnNeaNkvSfd31VMTAysVquydOjQAYAziPPz8zFv3jxMnDgR2dnZWL16NU6dOoW1a9fqXGoiIvJL09yVukdSSNhsNlRUVAAAVq9ejdOnT+tcovA5efIk+vTpg6VLl/rcHsjF8dmzZ2PDhg1Yt24dduzYgaqqKowdO5ZN44mIIpnWnDdI1nvSvU/1wYMHkZGRAYvFgoEDB2LBggU488wzUVhYiJKSEowcOVLZ12KxYMiQIdi5cydmzJihY6kp7Dyahnmtd9vcYNMwuc2Pqvm3Q3Ve4ev8nk3D5PO7zqc8N8d4NwWTp+rw9YVh0C+RYHm29nJv/tVQQzBHEC3ETG5Nveyuc3u2/pKbignBpmFh04JH/x40aBAmTJiA/v37QwiBO++8s8Hps15++eVmLl1o5eTkICcnx+c2z4vjgPMiQ1paGtauXYsZM2agvLwcL730El599VVcfvnlAIA1a9YgMzMTW7ZswahRo5rtvZAO/GW9xya/We/V/NtP1ntNfemd9aqcB5j1GrhnvWczb1+RrjXnAd9Z757zALM+LDj6t4qu72rgwIF45ZVX8NFHH+HFF19ESUkJBg8ejOPHjyv9qtPS0lTHuPe59qWmpgYVFRWqhYiImlkLvnq9Zs0ajB49GlVVVZAkCeXl5SgrK/O5GFljF8cBYPfu3aitrVXtk5GRgezsbGUfIiKKQLxTraLrnWr3q9u9e/fGoEGDcNZZZ2H16tW48MILAThHR3UnhPBa5y4vL89r8DMiIqLmkpaWhieffBIAkJWVhVdffRXt2rXTuVTNz9/F8UOHDin7xMXFoW3btl77NHYBvaamRnnOC+hERKSniLpU0KpVK/Tu3RsHDx6E1WoFAK9QLS0t9Qpod3PnzkV5ebmyHD58OKxlJiIib0KSNAxeYrz2eYWFhQFVqHv37m3YvAr24ngg++Tl5SElJUVZMjMzQ1JWIiIKjLacN2bWAxHQp9pdTU0NfvjhB1xyySXIysqC1WrF5s2b0a9fPwDOwV+2b9+OhQsXNngOi8UCi8XSXEWmZqJMreHZ78pk9u5vJfePcriuGTkc9f2f5J39TePgNh2W+2srUwC4vY5X/6oGyq2XSPjiauijlvtYCdU6+RjhsW8Qr+f2luU/yj37X8mvbZIkTrMVLi24T7UWv/76K2pra/UuRki5XxxPT09X1rtfHLdarbDZbCgrK1PdrS4tLcXgwYMbPPfcuXMxZ84c5XlFRQUr1gbgM+td+eo36+Wp60KV9Z6v4Sfr9c55Zxn0Dy5/We+5yVfWa815wHfWu+e8r/JFwEcW/dinWkXXd3Xvvfdi+/btKCwsxL///W9cc801qKiowJQpUyBJEmbPno0FCxZgw4YN+O677zB16lQkJiZi0qRJehabiIgaI0naFjIM94vjMvniuFxh7t+/P2JjY1X7FBcX47vvvvNbqbZYLEhOTlYtRETUjLTmvEGzXtc71UeOHMGNN96IY8eOoUOHDrjwwgvx5ZdfokuXLgCA+++/H9XV1Zg5cybKysowcOBAbNq0CUlJSXoWm4iIGsM71S1CVVUVfv75Z+V5YWEh9u3bh9TUVHTu3Fm5ON6tWzd069YNCxYsUF0cT0lJwfTp03HPPfegXbt2SE1Nxb333ovevXsro4ETEVEE4p1qFV0r1evWrfO7XZIk5ObmIjc3t3kKRJHHY7oNpWmY+z7yLkK9r6Q013LUT6fh2Xzc4dmc3O0/ekPTbrg/92g+JlM1BzPol0dj/E2dBdQ3+3aI+iZgSpMwt23Oc7k1EfN4Hc9PV5IkpQmYJDcxl+q3AfVNxNzLxOk3QkvuOxXsMRRddu3ahWHDhinP5SbZU6ZMwapVqwK6OP7ss88iJiYG1113HaqrqzF8+HCsWrUKZrPvbjVkQD6y3usr2F/WmzymzvKX9SaP7xkfWS+81jHrGxJI1ntmua+s9+r25fE67p+uZ5b7ynpfOe88hlkfKlpyXj7OiCKqTzURERFFj6FDh3r9MewukIvj8fHxWLJkCZYsWRKGEhIREYUfK9VERBR6ksn7jlAgxxAREVHk05Lz8nEGxEo1RYeGmoELh9t/TrmxkEczbeHwHjVcfu7r/3VDzbs8vwQkyX9TMF/HtBDuN648R/n21QzM7rFOfu6+r2dTME/yJy1J9U3CzB7Nvs1S403E2DQsRNinOij/+Mc//E4XSdQiuOW210jg/rJebgoe6qxns+8GNdTs21fWyz/7y3rlb4JGXteE+l+Lv6xvrNuXe9Yz5zVin2oVVqqJiCj0WKlWfPzxx/j4449RWloKh8c4Di+//DIAcFYLIiKKLqxUq7BSTUREocdKNQBg/vz5eOyxxzBgwACkp6crd02IiIiiGivVKqxUExFRyAlJ0jD6t/EqnMuXL8eqVaswefJkvYtCREQUMlpyXj7OiFippujSQN9qwG2aLfkHpR+WuX4KDk+e/a98fTl4/ud326fBL5MIuAqnx5dWMH2p5b5UdodQ/aza5tbvWj63A+qOXPLEK/V9rATMSr8q5zqz/M9Gee7aDgEH2N+Kwsdms2Hw4MF6F4Mourj1k/bMWZ9ZLzmnXwso6xvMbanBfSI16/WqnASb9V7Z7iPr649zncdP1ptc46P4y3r3nHeer+Gs5xgqFAr6/+VPRETGIzcLC3YxmL/85S9Yu3at3sUgIiIKLa05b8CsB1ipJiKicJAkbYsGL7zwArKyshAfH4/+/fvjs88+87v/9u3b0b9/f8THx+PMM8/E8uXLNb1uIE6fPo3FixdjyJAhuOOOOzBnzhzVQkREFJW05ryGrI/knJex+TdFJ/erXI01EXN/4tXc2xzYa7jx2QzMoFfdtHII4dUUTG7uJa+vlZt/OYA618/yujq3bYCzaVitXd0UzFOsWW4GJsFscu4b42rnFes6VH6uNCszSZBE/c9uD17NwClIzTRQ2fr16zF79my88MILuOiii/CPf/wDOTk5+P7779G5c2ev/QsLCzF69GjceuutWLNmDT7//HPMnDkTHTp0wNVXXx306zfmm2++Qd++fQEA3333nWobBy0j8sPz+6AZsz5Sm3tHkoaaffvKejnL/WW93N0rkKxXmn37yXr3nAfgN+uZ8xo100BlkZ7zMkkIz5nm/Js6dSqmTZuGSy+9NFxlCqmKigqkpKTg95ISJCcn610cCgfP8HTxmq/Sz76+TxDdlWo9+1TrX6l2rlOC1vUoP5e3x5ok5W8ws0fQyutNkmTYflYVFRWwpqWhvLw8ZN+Pyndu8W9Bn7OiogJp6RlBlWfgwIH405/+hGXLlinrzj33XEyYMAF5eXle+z/wwAN455138MMPPyjrbrvtNnz99df44osvgipvtIq2HA8Es74FaMasj5ZKtZ59qiOjUu1c5yvr3XMegN+slyvVzPrAz6c15+Xjg8n6aMn5oL8dKisrMXLkSHTr1g0LFizA0aNHw1EuIiKKZs3Qz8pms2H37t0YOXKkav3IkSOxc+dOn8d88cUXXvuPGjUKu3btQm1tbXDvMUoxx4mIqMmaoU91NOV80M2/33rrLRw/fhxr1qzBqlWr8Oijj+Lyyy/H9OnTMX78eMTGxoajnERE1EJUVFSonlssFlgsFq/9jh07BrvdjrS0NNX6tLQ0lJSU+Dx3SUmJz/3r6upw7NgxpKenN7H03goKCvDGG2+gqKgINptNte3tt98O+es1hjlORER6CyTroyXnAY0DlbVr1w533XUX9u7di6+++gpnn302Jk+ejIyMDNx99904ePBgqMtJ1LAGrnoJyeS1BHMVzdfxXk3CInAUQ+e8gc3bhkkI9RQbgLMpmHMqLOFc4GwOVusQqHUI1Nmdi83uwGnXUl3nXE7anMsf1bX4o7oWx07ZlKX0ZI1qcd927JQNf1TXKsfL55PPb3Mt8mvXOurLJZdTLndj74/8k/8dBrsAQGZmJlJSUpTFV/Mud559k4UQfvsr+9rf1/pQWLduHS666CJ8//332LBhA2pra/H999/jk08+QUpKSshfL1DMcYo6zZj1gb62nvTq4qUl6+Xs9Zf1nlnuL+vlvw38Zb17zjPrw0NrzmvJ+kjOeVmTBiorLi7Gpk2bsGnTJpjNZowePRr79+9Hz549sWjRItx9992hKicREUURLX+cyPsfPnxY1c/K111qAGjfvj3MZrPX1erS0lKvq9Qyq9Xqc/+YmBi0a9cuuAIHYMGCBXj22Wdx++23IykpCc899xyysrIwY8aMsF0tDwZznIiItNB6ESKYrI+GnJcFfdmttrYWb731FsaOHYsuXbrgjTfewN13343i4mKsXr0amzZtwquvvorHHnssHOUlIqIo4BBC0wIAycnJqqWhSnVcXBz69++PzZs3q9Zv3rwZgwcP9nnMoEGDvPbftGkTBgwYEJZmz//5z38wZswYAM4/GE6ePAlJknD33XdjxYoVIX+9QDDHiYioqbTmfDBZHw05Lwv6TnV6ejocDgduvPFGfPXVV8pUIe5GjRqFNm3ahKB4REFoqGmW2yigDY7o2ZTzt2DezcDqR/2Um9rIA3kqo326Vthcz212gdN1zt9RZY0dAHCqVv1YXWtXjq+1q0d1jXUN8SmP8JkQa0ZirHP6FPkxyeJ8FDEmpXxe5JHBXVvl5knu023I79eoI4SGktzULthjgjVnzhxMnjwZAwYMwKBBg7BixQoUFRXhtttuAwDMnTsXR48exSuvvALAOQLo0qVLMWfOHNx666344osv8NJLL+H111/X8OqNS01NRWVlJQDgjDPOwHfffYfevXvjxIkTOHXqVFheszHMcYpqvrLYY7RvZn1o+cp6oWxrPOttrp/9ZX216+dAsj7BI+N9ZX2DeeIj6z2n0GTWB0ZLzsvHBSPSc14WdKX62WefxbXXXov4+PgG92nbti0KCwubVDAiIopevvqrBXJMsK6//nocP34cjz32GIqLi5GdnY33338fXbp0AeBs3lxUVKTsn5WVhffffx933303nn/+eWRkZODvf/972OauvOSSS7B582b07t0b1113He666y588skn2Lx5M4YPHx6W12wMc5yIiJpKS87LxwUj0nNeFvQ81dGGc1dSUPNV+hMlV6+bc/ASf3eq5bkqI+lOdbzrTnWca07rOJOEGLN6fkvXU2UuS/c5LGVGuXodznmqi34L/ju3oqICnTOsIS2P3v744w+cPn0aGRkZcDgcePrpp7Fjxw6cffbZePjhh9G2bVu9i2gIzPoWLlQ5DzDrfb2WnzvVgWR9c9+pjo8xqXIegN+slz9JZn3g59Oa8/LxRst6oIkDlRFFhSgJyGjmXpl2Pq+/EqkErisfa90q0wBw0uZAla0OAFBe43o87XxU1lfXotpmV51PJld+E+KcYZqSEIvWcc6vtpT4GFX57K71iHP+m5BQPxKk3FhMkivTrpcxSfBqGkaNk0dYDfYYo0lNTVV+NplMuP/++3H//ffrWCIiA2LONwv3rJejOJCsP2lzrvSX9eXVzvmDA8n6lARnv1i/We+W84D/rHf9yKwPkpacl48zIlaqiYgo5Jqr+XckqqioUK6+e87D6clIV+mJiKjlaK7m39GClWoiIgoLg+Zmo9q2bYvi4mJ07NgRbdq08TkvpjzHpt1u16GERERETddSc94XVqqJiCjkWvKd6k8++URp9r1161adS0NERBR6vFOtFjGV6ry8PDz00EO46667kJ+fD8B5JX/+/PlYsWIFysrKMHDgQDz//PPo1auXvoUlikB6DlrivV0og5bIj3Vy/yqHesCSKlsd/nD1pTp+ygYAOHbS+fhHlfPxRHUtqk4796mpUw9eYnENPtY63tnHqk1CLFJbx7n2dT56981yfvVJkgkmSd2HSulL7WdODU630biW3Kd6yJAhPn8mImqqSMl6z6m0/GW93JfaX9afcG0LJOvbuPpU+8t695wH4D/rG/hMmfX+sU+1WkSM7FBQUIAVK1bgvPPOU61ftGgRFi9ejKVLl6KgoABWqxUjRoxQ5vwkIqLI5NC4GM2HH36IHTt2KM+ff/559O3bF5MmTUJZWZmOJWsedXV1+Nvf/oasrCwkJCTgzDPPxGOPPQaHo/63LYRAbm4uMjIykJCQgKFDh2L//v06lpqIiBqjNeeNmPVABFSqq6qqcNNNN+HFF19UTS0ihEB+fj7mzZuHiRMnIjs7G6tXr8apU6ewdu1aHUtMREQUmPvuu08ZrOzbb7/FnDlzMHr0aPzyyy+YM2eOzqULv4ULF2L58uVYunQpfvjhByxatAhPPfUUlixZouzDC+hERBTtdG/+ffvtt2PMmDG4/PLL8fjjjyvrCwsLUVJSgpEjRyrrLBYLhgwZgp07d2LGjBl6FJeIfHCfSkt+VJqGyXNOekylJc9TWV5TpzQF+72yBgBQfOK087G8GgBwosqGmmpn8zG7x9yVZtfclZYE59dZm9ZxOGVLUL22sq9r3gy5+ZfZJMHsagOmzFmpTK3hPd0GBU6IxrsJ+DrGaAoLC9GzZ08AwFtvvYVx48ZhwYIF2LNnD0aPHq1z6cLviy++wPjx4zFmzBgAQNeuXfH6669j165dALwvoAPA6tWrkZaWhrVr1zLriSKIr2kzA8l6eQotf1l/wtXdK5Csb+Nq9u0v691zHoDfrHefQpMCpyXn5eOMSNc71evWrcOePXuQl5fnta2kpAQAkJaWplqflpambPOlpqYGFRUVqoWIiJqX/EdXsIvRxMXF4dSpUwCALVu2KBeKU1NTW0Q+XXzxxfj444/x008/AQC+/vpr7NixQ7mg0NgFdCIiikxac96IWQ/oeKf68OHDuOuuu7Bp0ybEx8c3uJ/nVCTyNCQNycvLw/z580NWTiIiCl5LHqjM3cUXX4w5c+bgoosuwldffYX169cDAH766Sd06tRJ59KF3wMPPIDy8nL06NEDZrMZdrsdTzzxBG688UYA/i+gHzp0qMHz1tTUoKamRnneEi5QEBFFEg5Upqbbnerdu3ejtLQU/fv3R0xMDGJiYrB9+3b8/e9/R0xMjBKwnnelS0tLvcLX3dy5c1FeXq4shw8fDuv7ICIibxy8xGnp0qWIiYnBm2++iWXLluGMM84AAHzwwQe44oordC5d+K1fvx5r1qzB2rVrsWfPHqxevRpPP/00Vq9erdpPywX0lJQUZcnMzAxL+YmIyDcOVKam253q4cOH49tvv1Wtu+WWW9CjRw888MADOPPMM2G1WrF582b069cPAGCz2bB9+3YsXLiwwfNaLBZYLJawlp0okjTn9BqeHB5XG+WrjwJ+ptmocz6eqnX1qT5dp0yrIfevOnT8JACgrMz5vLqqJuA+1TXVdbB5TMUh96myxJidj65jEmPNiDW5+oO59pHLa3Z/T67PWH6/Js6v0SgBDX2qw1ISfXXu3Bn/+te/vNY/++yzOpSm+d1333148MEHccMNNwAAevfujUOHDiEvLw9TpkyB1WoF4LyAnp6erhwXyAV094HeKioqWLEmQ4u0rJfXBJL15aed+e0v66urnC1PAsl6eR9/We+e8wD8Zr3wmDuLWR8YLTkvH2dEulWqk5KSkJ2drVrXqlUrtGvXTlk/e/ZsLFiwAN26dUO3bt2wYMECJCYmYtKkSXoUmYiIKCh79uxBbGwsevfuDQD45z//iZUrV6Jnz57Izc1FXFycziUMr1OnTsFkUjeKM5vNypRaWVlZvIBORERRT/fRv/25//77UV1djZkzZ6KsrAwDBw7Epk2bkJSUpHfRiIjID4cQXnc3AjnGaGbMmIEHH3wQvXv3xi+//IIbbrgBV111Fd544w2cOnUK+fn5ehcxrMaNG4cnnngCnTt3Rq9evbB3714sXrwY06ZNA+Bs9s0L6ERE0UdLzsvHGVFEVaq3bdumei5JEnJzc5Gbm6tLeYhIG/fRHe3KdBvO53Kzb/mxylaHP1xTachTaMnNvqtcz0+W16DmZJXzPDXVqtcyW5xTalhatXZud2syFhfjav4V52z+lZIQCwBo7Xp+qtaOeNc+cjnlqbTk8pvZ+ksTgeCbeBkxZn/66Sf07dsXAPDGG2/g0ksvxdq1a/H555/jhhtuMHylesmSJXj44Ycxc+ZMlJaWIiMjAzNmzMAjjzyi7MML6ETRx30KTcB/1lfZnM21/WX9yXJX8+8Ast6zabivrHfPeQB+s545r42WnJePM6KIqlQTEZExaJk2w4jTbAghlKbOW7ZswdixYwEAmZmZOHbsmJ5FaxZJSUnIz8/3e/GAF9CJiKKP1umxjJj1ACvVREQUDkLDACYGDNoBAwbg8ccfx+WXX47t27dj2bJlAJzzM/sbiIuIiCiiacl513FGxEo1UZTSYyTQhr48PZuBuc9d6Nk0rNZ1165aHv27uhYnqmsBACdcTcPkUUDl5mCn/vgvaqr+cJ6ngSZh9ppU15oOyiihJ1wjgp9oXau8FgCkWGKUsng2BVNGMFfekwSpgabgHgOGkhsHBBxBJmew+0eD/Px83HTTTdi4cSPmzZuHs88+GwDw5ptvYvDgwTqXjogiXaRmfX1WNp71cvb6y/pTf/wXAALM+g7OdX6y3j3n3cvlK+vldcz64GjJefk4I2KlmoiIKEzOO+88r+kjAeCpp56C2Wz2cQQRERFFG1aqiYgo5ISGZmEGHRDUp/j4eL2LQEREpJmWnJePMyJWqomIKOQ4UJmTyWSC5KfNoN1ub8bSEBERhQYHKlNjpZqIwkKe8EK+IllrF+pH17dqtc2OqtPOflA11XXqR9fUGjVVf6Cm3DlScp1HP6sYVz8rmdmSgBrX1FnyeeTzV9vsqteutQulPCJGXW5qGt6pdtqwYYPqeW1tLfbu3YvVq1dj/vz5OpWKiCg0Asl6OXv9Zb3clzqQrJf7V/vLevecd39k1ocO71SrsVJNREQhx4HKnMaPH++17pprrkGvXr2wfv16TJ8+XYdSERERNQ0HKlMz6V0AIiIyHvkKdrBLSzFw4EBs2bJF72IQERFpojXnjZr1vFNNRCHnUP2s/vZ0yNNt2F3TXDgEaupcP9s9Hl3Nv+w11UpTMLtN3SRMFuO2r92eojqPcn6H+rUdbt/s9eWUvN4DUShVV1djyZIl6NSpk95FISLSxDMj/WW9nL3+sl7O+0CyXvnbwE/W+8p5dTmZ9RRarFQTEVHIOYTw+mMmkGOMpm3btqqByoQQqKysRGJiItasWaNjyYiIiLTTkvPycUbESjUREYWc3eFcgj3GaPLz81XPTSYTOnTogIEDB6Jt27b6FIqIiKiJtOS8fJwRsVJNREQhxzvVTlOmTNG7CERERCHHO9VqHKiMKMoISYLwM+8tUSRwCAF7kEs4g7asrAyTJ09GSkoKUlJSMHnyZJw4ccLvMVOnToUkSarlwgsv1PT6p06dwo8//ohvvvlGtRAR+cKsp0inJeeNnPW8U01ERCHnEMFfjXaE8eL1pEmTcOTIEXz44YcAgP/3//4fJk+ejHfffdfvcVdccQVWrlypPI+Liwvqdf/73//illtuwQcffOBzu91uD+p8REREkUBLzsvHhYteWQ+wUk1ERAb3ww8/4MMPP8SXX36JgQMHAgBefPFFDBo0CAcOHMA555zT4LEWiwVWq1Xza8+ePRtlZWX48ssvMWzYMGzYsAG///47Hn/8cTzzzDOaz0tERET19Mx6gJVqIgoDk+pndfM1k6s5W6zZuZfZJMES4/rZ7PFoSVAeY1w/e4px20d+9DyPcn6T+rVNbk3rvMrZ+NskP5oyUFlFRYVqvcVigcVi0VyWL774AikpKUrIAsCFF16IlJQU7Ny502/Qbtu2DR07dkSbNm0wZMgQPPHEE+jYsWPAr/3JJ5/gn//8J84//3yYTCZ06dIFI0aMQHJyMvLy8jBmzBjN74uISC+eGekv6+Xs9Zf1coYHkvVK3vvJel8577Ocft8l+dPUgcqMlPUA/y0REVEYyAOYBLsAQGZmptIfKiUlBXl5eU0qS0lJic9w7NixI0pKSho8LicnB6+99ho++eQTPPPMMygoKMBll12GmpqagF/75MmTymunpqbiv//9LwCgd+/e2LNnT5DvhIiIKDJozXkjZj3AO9VERBQG8oAkwR4DAIcPH0ZycrKyvqEr17m5uZg/f77fcxYUFACAaq5omRDC53rZ9ddfr/ycnZ2NAQMGoEuXLnjvvfcwceJEv68rO+ecc3DgwAF07doVffv2xT/+8Q907doVy5cvR3p6ekDnICIiijRacl4+DjBW1gOsVBNFlWgaCVRuBiMXOdYsqR9dzcES4sxoHR8LALAkxKgfW7UGANhrUpXzxtRUq15HbgZmaZ2qHON5Hvn8CXFm1WvHmiWlPHI52XwnNBwIfjASuRVZcnKyKmgbMmvWLNxwww1+9+natSu++eYb/P77717b/vvf/yItLS3g8qWnp6NLly44ePBgwMfMnj0bxcXFAIBHH30Uo0aNwpo1axAXF4fVq1cHfB4iahmiKeeBwLJezl5/We+e84D/rJf/NvCX9e457/7IrA8dLTkvHwcYK+sBVqqJiCgM7A4Be5BpG+z+7du3R/v27Rvdb9CgQSgvL8dXX32FCy64AADw73//G+Xl5Rg8eHDAr3f8+HEcPnw4qDvMN910k/Jzv3798Ouvv+LHH39E586dAyo7ERFRJNKS8/JxwYiGrAdYqSYiIoM799xzccUVV+DWW2/FP/7xDwDOaTbGjh2rGrikR48eyMvLw1VXXYWqqirk5ubi6quvRnp6On799Vc89NBDaN++Pa666qqAX3vOnDk+10uShPj4eJx99tkYP348UlNTfe5HREREjdMz6wFWqomIKAyE22AkwRwTLq+99hruvPNOjBw5EgBw5ZVXYunSpap9Dhw4gPLycgCA2WzGt99+i1deeQUnTpxAeno6hg0bhvXr1yMpKSng1927dy/27NkDu92Oc845B0IIHDx4EGazGT169MALL7yAe+65Bzt27EDPnj1D94aJiIjCSEvOy8eFi15ZD+hcqV62bBmWLVuGX3/9FQDQq1cvPPLII8jJyQHg/NDnz5+PFStWoKysDAMHDsTzzz+PXr166VhqopZL7ovk+X1oktdL8n5Sfb8lybmzWZ5ew+TsyZQQ6+xjlZIQizYJzn5QbVrHAQBqqusAAHZlroYOSn8qe0P9rFx9rFqlWJDQ2qI6n3z+FNej/NqxJpNSLpNb2d0fTVL9toY+D/JmF84l2GPCJTU1FWvWrPG7j3vQJyQk4KOPPmry68p3oVeuXKn0HauoqMD06dNx8cUX49Zbb8WkSZNw9913h+T1iIiaKpCsr89KeVvDWS9nr/+s7+A8PoCsb5XizHh/We+e8+7l8pX1JuU9+P88SE1LzsvHhYteWQ/o3E+/U6dOePLJJ7Fr1y7s2rULl112GcaPH4/9+/cDABYtWoTFixdj6dKlKCgogNVqxYgRI1BZWalnsYmIqBFNmWbDSJ566in8z//8j2owluTkZOTm5mLRokVITEzEI488gt27d+tYyvA6evQo/vznP6Ndu3ZITExE3759Ve9XCIHc3FxkZGQgISEBQ4cOVf4OICKiyNTUKbWMRtdK9bhx4zB69Gh0794d3bt3xxNPPIHWrVvjyy+/hBAC+fn5mDdvHiZOnIjs7GysXr0ap06dwtq1a/UsNhERNUIewCTYxWjKy8tRWlrqtf6///0vKioqAABt2rSBzWZr7qI1i7KyMlx00UWIjY3FBx98gO+//x7PPPMM2rRpo+zDC+hERNFHa84bMeuBCOpTbbfb8cYbb+DkyZMYNGgQCgsLUVJSorSJB5zzlw0ZMgQ7d+7EjBkzdCwtEfnj3mRabnJlNjm/RBNdTbLkx9ZxMUh1Nd06ZXM277LVOdxPB7PZhBpXsy67PcVrG1A/pUZCawvato0HAKSnOM8nn791XIxXGVyHezUNa6gZGAVGy9VoI169Hj9+PKZNm4ZnnnkG559/PiRJwldffYV7770XEyZMAAB89dVX6N69u74FDZOFCxciMzMTK1euVNZ17dpV+dnzAjoArF69GmlpaVi7di2znihCeWalv6yXs9df1stZHkjWy82+/WW9598azPrQ03rX2YhZD0TANG3ffvstWrduDYvFgttuuw0bNmxAz549UVJSAgBe84qlpaUp23ypqalBRUWFaiEiItLDP/7xDwwfPhw33HADunTpgs6dO+OGG27A8OHDsXz5cgDOkUj/93//V+eShsc777yDAQMG4Nprr0XHjh3Rr18/vPjii8r2xi6gExERRQPd71Sfc8452LdvH06cOIG33noLU6ZMwfbt25XtksfoAEIIr3Xu8vLyMH/+/LCVl4iIGhdpA5XppXXr1njxxRfx7LPP4pdffoEQAmeddRZat26t7NO3b1/9Chhmv/zyC5YtW4Y5c+bgoYcewldffYU777wTFosFN998s98L6IcOHWrwvDU1NaipqVGe8wI6EVHzisSByvSke6U6Li4OZ599NgBgwIABKCgowHPPPYcHHngAAFBSUqKafLu0tNQrfN3NnTtXNS9oRUUFMjMzw1R6IiLyhc2/1Vq3bo3zzjtP72I0O4fDgQEDBmDBggUAgH79+mH//v1YtmwZbr75ZmU/XkAnIooubP6tpnul2pMQAjU1NcjKyoLVasXmzZvRr18/AIDNZsP27duxcOHCBo+3WCywWCzNVVyiZiEidD4Hk6tcdtcXpDJFhRAwK/2snI8xro5LcTHOR7mfU0p8DGrqnP2gPAeviItx9lA5kRDjY5ot1/k9+lS3aR2n9K9Kb+Psb9W+VZzyWu6vHRcjKeXyLK/8ibv/YW+K0N9DJHI4BBxBDkYS7P4U+dLT073m3z733HPx1ltvAQCsVisAXkAnAqIr6yUhT6EF1aOvrJez11/Wn3BleCBZL0+h5S/r3XPevVy+st7zAh6zPjBacl4+zoh0rVQ/9NBDyMnJQWZmJiorK7Fu3Tps27YNH374ISRJwuzZs7FgwQJ069YN3bp1w4IFC5CYmIhJkybpWWwiImqEQ0OzMIPmbIt20UUX4cCBA6p1P/30E7p06QIAvIBORBSltOS8fJwR6Vqp/v333zF58mQUFxcjJSUF5513Hj788EOMGDECAHD//fejuroaM2fORFlZGQYOHIhNmzYhKSlJz2ITEVEj2PybAODuu+/G4MGDsWDBAlx33XX46quvsGLFCqxYsQIAeAGdiChKsfm3mq6V6pdeesnvdkmSkJubi9zc3OYpEBFpIjeUkqemEFJ9cyp5GotY13eo3dXmKsnibJrlEN5zFppdJ0qMc+5zonUtqk7XAgBqPKbbsriaiLeOd07D0SYhVplWQ24K1i7R1STMEqN67TizhFi5SZjrUfIx3QYbghFpc/7552PDhg2YO3cuHnvsMWRlZSE/Px833XSTsg8voBNFB/esF3Lz6QCyXq5E+cv6E62dGR9I1rdxTbvlL+vdcx6A36w3uTUFJ9Iq4vpUExFR9LMLofS/C+YYMp6xY8di7NixDW7nBXQiouijJefl44yIlWoiIgo5DlRGRERkXByoTI2VaqIIFqkjgcrFauhioyRJMEvOjXITMWVEUNcK4WrKZY+r/xqSm2VZYlwjg7uaeJVX16LaZnfu30DzsYS4+mNau84pjzgqN/uW18e7XjvO5Gf0bz+ffYT+WiKKHRrmqQ5LSYiIIls0Zr3S7DuArHfPecB31pdXO5t9B5L18t8G/rLePefdy8WsDx0tOS8fZ0SsVBMRUchxoDIiIiLj4kBlaqxUExFRyLFPNRERkXGxT7WaSe8CEBEREREREUUr3qkmoiYzuToeyU163KfbkPtBOeDcFiv3s5IPjjPBbIpRncfimpujtdxP2hKDWlf/qlq7epqNWNe+8nQZCbFmJMY6j5Mf5ak1lD5WblNsyNOAyOV0n0oLcE6xYWLHqqA5HN5TpQVyDBERRSb3rHefQhNoJOvjTK59Gs56uT90IFmf4JHxvrLeeyotdTl9TZvJrA+OlpyXjzMiVqqJiCjk7BrCVks4ExERUfPTkvPycUbESjUREYUcK9VERETGxUq1GivVRBEoUqfX8OQ53YZ70zBlug0oG32fA4AkyU3DnPvITblO1TonXqh1OFDrmrfBc9RI+TVjlaZeJuX4uBjX1B6ubfLUGnJzsBizVN80zGN6DV/NwaLk1xIR7I7gg9OjtR8RkaFFc9bLWRxI1str/GV9rcMZAIFkfazJpDreV9a75zwAv1nv2ew7Sn4tutOS8/JxRsRKNRERhRzvVBMRERkX71SrcfRvIiIiIiIiIo14p5qIiEKOd6qJiIiMi3eq1VipJqIm8+xvBdR3q3K4elNJro2xJs++SxJMknObWZL3cT7KU2DZhVD6WTVE7lNtluqnyYrx7EPtMbVGrEly6+ulnmbD1/ujwHFKLSIiY9Ga9fX52nDW24U8lVbjWW+W1FnuK+vdcx4Asz4MOKWWGivVREQUcnah4U61MGbQEhERGY2WnJePMyJWqomIKOTY/JuIiMi42PxbjZVqIgo59+k2lCZWrh9E/ZwczgcIZToLeaoLu2ubfDXTISQI17eVA+ovY5Pc5EyqP63Z43xmU8PPPZuC+ZpKi4LHSjURkXG5T6HpfK5sAKDOehOE6hhfWe8QruMCyHr5tfxlvWfuM+tDj5VqNY7+TUREhvfEE09g8ODBSExMRJs2bQI6RgiB3NxcZGRkICEhAUOHDsX+/fvDW1AiIiLSRM+sZ6WaiIhCrs4hNC3hYrPZcO211+Kvf/1rwMcsWrQIixcvxtKlS1FQUACr1YoRI0agsrIybOUkIiKKBlpz3qhZz+bfRBFERHkzJPeRQRtqGiaPEGp2Ne2STBJcA4HWNwlzPVeagwkBR/2r+Hxt+QqhczRx9fkkjyZiJsl734aagkX5r0Q3kdb8e/78+QCAVatWBbS/EAL5+fmYN28eJk6cCABYvXo10tLSsHbtWsyYMSNcRSUigzNC1te37m486yV5BG4/WS83Fw8k6z2bcvvKepPHOn9ZH+W/Dt1EYvNvPbOed6qJiCjk5Kk2glkiaZqNwsJClJSUYOTIkco6i8WCIUOGYOfOnTqWjIiISH9act7IWc871UREFHJ2IYKeNkPev6KiQrXeYrHAYrGErGyBKCkpAQCkpaWp1qelpeHQoUPNWhYiIqJIoyXn5eMA42U971QTEVHIabl6LTcJy8zMREpKirLk5eX5fI3c3FxIkuR32bVrV5Peh+TRLlAI4bWOiIiopdGa80bNel3vVOfl5eHtt9/Gjz/+iISEBAwePBgLFy7EOeeco+wjhMD8+fOxYsUKlJWVYeDAgXj++efRq1cvHUtORP74628l932S+5RJAjC5TZ0FAGbXeeQWQu4XQuv7W7nO7/XavvpJ129TPXc/D/tSR4zDhw8jOTlZed7QletZs2bhhhtu8Huurl27aiqD1WoF4LyKnZ6erqwvLS31uqJNRNTSuI+hAvjPerkvtb+s97zh6S/rG+4nXb/dV867l9PzfVDzM1rW63qnevv27bj99tvx5ZdfYvPmzairq8PIkSNx8uRJZR+OvkpEFH2acvU6OTlZtTQUtO3bt0ePHj38LvHx8ZrKn5WVBavVis2bNyvrbDYbtm/fjsGDB2s6JzkvpkuShNmzZyvrOHUZEVH0aeqdaqNlva6V6g8//BBTp05Fr1690KdPH6xcuRJFRUXYvXs3AO8R2bKzs7F69WqcOnUKa9eu1bPoRETkR1OCNhyKioqwb98+FBUVwW63Y9++fdi3bx+qqqqUfXr06IENGzYAgFLxW7BgATZs2IDvvvsOU6dORWJiIiZNmhS2chpZQUEBVqxYgfPOO0+1nhfPiYiiT1Mr1eGgZ9ZH1EBl5eXlAIDU1FQAjY/IxilNyEiifYoNTw01DZOppt9wawruPEbdfExuKhYIk9uunk3ElPXKvt7nNdivQTd24YDd4dmAr/FjwuWRRx7B6tWrlef9+vUDAGzduhVDhw4FABw4cEDJIQC4//77UV1djZkzZyrdjzZt2oSkpKSwldOoqqqqcNNNN+HFF1/E448/rqzn1GXUkhgt54HAsl7JYD9ZrzXnnaf1znp/Oe9ebtJOS87Lx4WLnlkfMZVqIQTmzJmDiy++GNnZ2QC0jchWU1ODmpoa5bnnyHJERBR+Dg1Xo8M5zcaqVasanbdSeHTqkyQJubm5yM3NDVu5Worbb78dY8aMweWXX66qVPPiORFRdNKS8/Jx4aJn1kdMpXrWrFn45ptvsGPHDq9twYzIlpeXp0z8TURE+rA7BExBBmc4m4SRftatW4c9e/agoKDAa5vW6Ux4AZ2ISF9acl4+zogiYkqtO+64A++88w62bt2KTp06KevdR2Rz529Etrlz56K8vFxZDh8+HL6CExERUYMOHz6Mu+66C2vWrPE7kEyw05nk5eWppmLJzMwMWZmJiIiCpWulWgiBWbNm4e2338Ynn3yCrKws1XYtI7JZLBav0eSISD+SpF5kJkmCSZIgAcpikpyL2SSFZJHP5/4akttr+ysfNU2dA6hziCAXvUtNobZ7926Ulpaif//+iImJQUxMDLZv346///3viImJUS6QB3PxHOAFdKJI4y/rvTM4/FnvnvP+ykfaact542a9rs2/b7/9dqxduxb//Oc/kZSUpIRqSkoKEhISVCOydevWDd26dcOCBQs4+ioRUYRj828CgOHDh+Pbb79VrbvlllvQo0cPPPDAAzjzzDOVi+fygDLyxfOFCxc2eF6LxdLg9CtERBR+bP6tpmuletmyZQCgjMYmW7lyJaZOnQqAo68SEUUjVqoJAJKSkpTBR2WtWrVCu3btlPW8eE5EFH1YqVbTtVLtOfqaLxx9lchYGpt+A6ifbktm1thUq6GpNDzLQqHHSjUFihfPiYzHPesbymL3rNea84D/rGfOhw8r1WoRM/o3EREZR6RNqUWRY9u2barnvHhORBR9InFKLT2xUk1ERCFndwhIvFNNRERkSFpyXj7OiFipJtKZaKFtk9zftmdPkMaabYfiNYmIiJoLs55ZT8bGSjUREYWcEAIiyKvRgYyzQURERPrTkvPycUbESjUREYWcwyGC7jdl1H5WRERERqMl5+XjjIiVaiIiCjkhRNBXo4169ZqIiMhotOS8fJwRsVJNRLprrP9TMN+/7EsVGYRDQ/Nvg169JiIi//nMnI8+WnJePs6IWKkmIqKQY/NvIiIi42LzbzWT3gUgIiIiIiIiila8U02kk5Y6vYYW/Kiij3A4l2CPISIyEmZ9YPgxRR8tOS8fZ0SsVBMRUchxoDIiIiLj4kBlaqxUExFRyLFPNRERkXGxT7UaK9VERBRyHP2biIjIuDj6txor1UREFHpawtagQUtERGQ4GivVRs16jv5NREREREREpBHvVBMRUcg5hIAU5GAkDoMOXkJERGQ0WnJePs6IWKkmamacXoNaAiE09Kk2aNASUcvDrCej05Lz8nFGxEo1ERGFHAcqIyIiMi4OVKbGSjUREYWcwwFIQU+pFabCEBERUUhpyXn5OCNipZqIiEJOCBF0Ey+jNgkjIiIyGi05Lx9nRBz9m4iIiIiIiEgj3qkmIqKQEw7nEuwxREREFPm05Lx8nBHxTjUREYWcwyE0LeHyxBNPYPDgwUhMTESbNm0COmbq1KmQJEm1XHjhhWErIxERUbTQmvNGzXpdK9Wffvopxo0bh4yMDEiShI0bN6q2CyGQm5uLjIwMJCQkYOjQodi/f78+hSUiooDJo4IGu4SLzWbDtddei7/+9a9BHXfFFVeguLhYWd5///0wlZCIiCh6aM15o2a9rpXqkydPok+fPli6dKnP7YsWLcLixYuxdOlSFBQUwGq1YsSIEaisrGzmkhIRUTAiLWjnz5+Pu+++G7179w7qOIvFAqvVqiypqalhKqEx5eXl4fzzz0dSUhI6duyICRMm4MCBA6p9eAGdiCj6RGKlWs+s17VSnZOTg8cffxwTJ0702iaEQH5+PubNm4eJEyciOzsbq1evxqlTp7B27VodSktERIFyCKFpAYCKigrVUlNTo9v72LZtGzp27Iju3bvj1ltvRWlpqW5liUbbt2/H7bffji+//BKbN29GXV0dRo4ciZMnTyr78AI6EVH00ZrzRs36iO1TXVhYiJKSEowcOVJZZ7FYMGTIEOzcuVPHkhERUThlZmYiJSVFWfLy8nQpR05ODl577TV88skneOaZZ1BQUIDLLrtM1+CPNh9++CGmTp2KXr16oU+fPli5ciWKioqwe/duALyATkTUUhkt6yN29O+SkhIAQFpammp9WloaDh061OBxNTU1qg+hoqIiPAUkIqIGaWniJe9/+PBhJCcnK+stFovP/XNzczF//ny/5ywoKMCAAQOCKofs+uuvV37Ozs7GgAED0KVLF7z33ns+W1hR48rLywFAaVrX2AX0GTNm6FJOIiLyT2tTbqNmfcRWqmWSJKmeCyG81rnLy8tr9IMnIqLwEkJDpdrVJCw5OVkVtA2ZNWsWbrjhBr/7dO3aNagy+JOeno4uXbrg4MGDITtnSyKEwJw5c3DxxRcjOzsbAC+gExFFKy05Lx8HGC/rI7ZSbbVaATgDNz09XVlfWlrqFb7u5s6dizlz5ijPKyoqkJmZGb6CEhGRF6Fh2oxgw7l9+/Zo3759UMc0xfHjx3H48GFVJlHgZs2ahW+++QY7duzw2sYL6ERE0UVLzsvHBSNasj5i+1RnZWXBarVi8+bNyjqbzYbt27dj8ODBDR5nsViUKx+BXgEhIqLQEkJoWsKlqKgI+/btQ1FREex2O/bt24d9+/ahqqpK2adHjx7YsGEDAKCqqgr33nsvvvjiC/z666/Ytm0bxo0bh/bt2+Oqq64KWzmN6o477sA777yDrVu3olOnTsp69wvo7gK5gF5eXq4shw8fDk/BiYjIJ605b9Ss1/VOdVVVFX7++WfleWFhIfbt24fU1FR07twZs2fPxoIFC9CtWzd069YNCxYsQGJiIiZNmqRjqYmIqDFN6VMdDo888ghWr16tPO/Xrx8AYOvWrRg6dCgA4MCBA0qfX7PZjG+//RavvPIKTpw4gfT0dAwbNgzr169HUlJS2MppNEII3HHHHdiwYQO2bduGrKws1Xb3C+jy70S+gL5w4cIGz2uxWBrsf0dEROHX1D7V4aBn1utaqd61axeGDRumPJebbU+ZMgWrVq3C/fffj+rqasycORNlZWUYOHAgNm3axD9oiIgoKKtWrcKqVav87uN+9TwhIQEfffRRmEtlfLfffjvWrl2Lf/7zn0hKSlLuSKekpCAhIQGSJPECOhERhYSeWa9rpXro0KF+mwBIkoTc3Fzk5uY2X6GIiKjJHA4BBHk1WkvfLIpsy5YtAwDlDoFs5cqVmDp1KgDwAjoRURTSkvPKcQYUsQOVERFR9BIOO4TDHvQxZCyB9J3jBXQiouijJefl44yIlWoiIgo5VqqJiIiMi5VqNVaqiZqZ5LpzI/xMF0MUbpIQyr/FcBAOh4ZKtSNMpSEial7MeooE4cx6LTkvH2dErFQTEVHICbsdwh5kpTrI/YmIiEgfWnJePs6IInaeaiIiIiIiIqJIxzvVREQUckJo6FMtjHn1moiIyGi05Lx8nBGxUk2kE/a3Ir2Esy+1jAOVEREx60kfkZrz8nFGxEo1ERGFHCvVRERExsVKtRor1UREFHKsVBMRERkXK9VqrFQTEVHIcUotIiIi4+KUWmoc/ZuIiIiIiIhII96pJiKikHM47ECQV7AdBm0SRkREZDRacl45zoBYqSYiopBjn2oiIiLjYp9qNVaqiXTG6TaouTTHFBsyVqqJiOpJQjDnKewiPefl44yIlWoiIgo9ux3CFGRw2o0ZtERERIajJeddxxkRK9VERBRyQgTf10oIYwYtERGR0WjJeeU4A+Lo30QRojmb7FDLIgnBf19ERDrjdzGFE/9t6Yt3qomIKOSEwxH8nWqDzl1JRERkNFpyXjnOgFipJiKikBMaptow6uAlRERERqMl55XjDIiVaiIiCjnnFezgrkYb9eo1ERGR0WjJeeU4A2KlmiiCcHotCiU9+1fxTjURkW/MegolvbKed6rVWKkmIqKQY6WaiIjIuFipVuPo30REREREREQa8U41ERGFnMNhh8Q71URERIakJecB42Z9VNypfuGFF5CVlYX4+Hj0798fn332md5FIiIiP4TdAWG3B7kYc/ASCgyznogoemjLeeNmfcRXqtevX4/Zs2dj3rx52Lt3Ly655BLk5OSgqKhI76IREVEDhLBDOIJchDGvXlPjmPVERNFFU84bOOsjvlK9ePFiTJ8+HX/5y19w7rnnIj8/H5mZmVi2bJneRSMKG0kIXUdupugWCf9+NAWtQZuEUeOY9dQSRcJ3NUUvvf/9aM15o2Z9RFeqbTYbdu/ejZEjR6rWjxw5Ejt37vR5TE1NDSoqKlQLERE1LwYtBUpL1hMRkb5YqVaL6Er1sWPHYLfbkZaWplqflpaGkpISn8fk5eUhJSVFWTIzM5ujqERERKSBlqznBXQiIookUTH6tyRJqudCCK91srlz52LOnDnK8/LycnTu3BmVlZVhLSNROIgG/p0T+RNoczD5e1GEofmYqD0d/NVoe23Iy0HRI5isz8vLw/z5873WM+spGjHrSQu9s15TzgOGzfqIrlS3b98eZrPZ60p1aWmp1xVtmcVigcViUZ7LV6/P7tYtfAUlIopilZWVSElJCcm54uLiYLVaUfL9/2k63mq1Ii4uLiRloeigJes9L6AXFhaib9++zHoiogaEKuubmvOAMbM+oivVcXFx6N+/PzZv3oyrrrpKWb9582aMHz8+oHNkZGTg+++/R8+ePXH48GEkJyeHq7gE50WMzMxMftbNgJ918zHqZy2EQGVlJTIyMkJ2zvj4eBQWFsJms2k6Pi4uDvHx8SErD0U+LVnveQG9S5cuAICioqKQXSAi34z6fRip+Hk3H6N+1qHO+qbmPGDMrI/oSjUAzJkzB5MnT8aAAQMwaNAgrFixAkVFRbjtttsCOt5kMuGMM84AACQnJxvqP0kk42fdfPhZNx8jftbhqIDEx8cbLiwpvEKR9YDz37PR/o9GKiN+H0Yyft7Nx4ifdaiznjnvLeIr1ddffz2OHz+Oxx57DMXFxcjOzsb777+vXJUmIiKi6MasJyKiaBbxlWoAmDlzJmbOnKl3MYiIiChMmPVERBStInpKrVCxWCx49NFHVf2vKDz4WTcfftbNh581UWTj/9Hmw8+6efHzbj78rKkpJBGOuVSIiIiIiIiIWoAWcaeaiIiIiIiIKBxYqSYiIiIiIiLSiJVqIiIiIiIiIo0MX6l+4YUXkJWVhfj4ePTv3x+fffaZ3kWKerm5uZAkSbVYrVZluxACubm5yMjIQEJCAoYOHYr9+/frWOLo8emnn2LcuHHIyMiAJEnYuHGjansgn21NTQ3uuOMOtG/fHq1atcKVV16JI0eONOO7iB6Nfd5Tp071+rd+4YUXqvbh502kP2Z90zF/mkdeXh7OP/98JCUloWPHjpgwYQIOHDig2oefdWgsW7YM5513njLv9KBBg/DBBx8o2/k5UygZulK9fv16zJ49G/PmzcPevXtxySWXICcnB0VFRXoXLer16tULxcXFyvLtt98q2xYtWoTFixdj6dKlKCgogNVqxYgRI1BZWaljiaPDyZMn0adPHyxdutTn9kA+29mzZ2PDhg1Yt24dduzYgaqqKowdOxZ2u7253kbUaOzzBoArrrhC9W/9/fffV23n502kL2Z9aDB/msf27dtx++2348svv8TmzZtRV1eHkSNH4uTJk8o+/KxDo1OnTnjyySexa9cu7Nq1C5dddhnGjx+vVJz5OVNICQO74IILxG233aZa16NHD/Hggw/qVCJjePTRR0WfPn18bnM4HMJqtYonn3xSWXf69GmRkpIili9f3kwlNAYAYsOGDcrzQD7bEydOiNjYWLFu3Tpln6NHjwqTySQ+/PDDZit7NPL8vIUQYsqUKWL8+PENHsPPm0h/zPrQY/40n9LSUgFAbN++XQjBzzrc2rZtK/73f/+XnzOFnGHvVNtsNuzevRsjR45UrR85ciR27typU6mM4+DBg8jIyEBWVhZuuOEG/PLLLwCAwsJClJSUqD53i8WCIUOG8HNvokA+2927d6O2tla1T0ZGBrKzs/n5a7Rt2zZ07NgR3bt3x6233orS0lJlGz9vIn0x65sH8yd8ysvLAQCpqakA+FmHi91ux7p163Dy5EkMGjSInzOFnGEr1ceOHYPdbkdaWppqfVpaGkpKSnQqlTEMHDgQr7zyCj766CO8+OKLKCkpweDBg3H8+HHls+XnHnqBfLYlJSWIi4tD27ZtG9yHApeTk4PXXnsNn3zyCZ555hkUFBTgsssuQ01NDQB+3kR6Y9Y3D+ZPeAghMGfOHFx88cXIzs4GwM861L799lu0bt0aFosFt912GzZs2ICePXvyc6aQi9G7AOEmSZLquRDCax0FJycnR/m5d+/eGDRoEM466yysXr1aGcSJn3v4aPls+flrc/311ys/Z2dnY8CAAejSpQvee+89TJw4scHj+HkTNS9mTvNg/oTWrFmz8M0332DHjh1e2/hZh8Y555yDffv24cSJE3jrrbcwZcoUbN++XdnOz5lCxbB3qtu3bw+z2ex1Jam0tNTrqhQ1TatWrdC7d28cPHhQGQWcn3voBfLZWq1W2Gw2lJWVNbgPaZeeno4uXbrg4MGDAPh5E+mNWd88mD+hd8cdd+Cdd97B1q1b0alTJ2U9P+vQiouLw9lnn40BAwYgLy8Pffr0wXPPPcfPmULOsJXquLg49O/fH5s3b1at37x5MwYPHqxTqYyppqYGP/zwA9LT05GVlQWr1ar63G02G7Zv387PvYkC+Wz79++P2NhY1T7FxcX47rvv+PmHwPHjx3H48GGkp6cD4OdNpDdmffNg/oSOEAKzZs3C22+/jU8++QRZWVmq7fysw0sIgZqaGn7OFHo6DI7WbNatWydiY2PFSy+9JL7//nsxe/Zs0apVK/Hrr7/qXbSods8994ht27aJX375RXz55Zdi7NixIikpSflcn3zySZGSkiLefvtt8e2334obb7xRpKeni4qKCp1LHvkqKyvF3r17xd69ewUAsXjxYrF3715x6NAhIURgn+1tt90mOnXqJLZs2SL27NkjLrvsMtGnTx9RV1en19uKWP4+78rKSnHPPfeInTt3isLCQrF161YxaNAgccYZZ/DzJoogzPrQYP40j7/+9a8iJSVFbNu2TRQXFyvLqVOnlH34WYfG3LlzxaeffioKCwvFN998Ix566CFhMpnEpk2bhBD8nCm0DF2pFkKI559/XnTp0kXExcWJP/3pT8qUBaTd9ddfL9LT00VsbKzIyMgQEydOFPv371e2OxwO8eijjwqr1SosFou49NJLxbfffqtjiaPH1q1bBQCvZcqUKUKIwD7b6upqMWvWLJGamioSEhLE2LFjRVFRkQ7vJvL5+7xPnTolRo4cKTp06CBiY2NF586dxZQpU7w+S37eRPpj1jcd86d5+PqMAYiVK1cq+/CzDo1p06Yp3wsdOnQQw4cPVyrUQvBzptCShBCi+e6LExERERERERmHYftUExEREREREYUbK9VEREREREREGrFSTURERERERKQRK9VEREREREREGrFSTURERERERKQRK9VEREREREREGrFSTURERERERKQRK9VEREREREREGrFSTURERERERKQRK9VEREREREREGrFSTURERERERKQRK9VEYfbf//4XVqsVCxYsUNb9+9//RlxcHDZt2qRjyYiIiCgUmPVELZskhBB6F4LI6N5//31MmDABO3fuRI8ePdCvXz+MGTMG+fn5eheNiIiIQoBZT9RysVJN1Exuv/12bNmyBeeffz6+/vprFBQUID4+Xu9iERERUYgw64laJlaqiZpJdXU1srOzcfjwYezatQvnnXee3kUiIiKiEGLWE7VM7FNN1Ex++eUX/Pbbb3A4HDh06JDexSEiIqIQY9YTtUy8U03UDGw2Gy644AL07dsXPXr0wOLFi/Htt98iLS1N76IRERFRCDDriVouVqqJmsF9992HN998E19//TVat26NYcOGISkpCf/617/0LhoRERGFALOeqOVi82+iMNu2bRvy8/Px6quvIjk5GSaTCa+++ip27NiBZcuW6V08IiIiaiJmPVHLxjvVRERERERERBrxTjURERERERGRRqxUExEREREREWnESjURERERERGRRqxUExEREREREWnESjURERERERGRRqxUExEREREREWnESjURERERERGRRqxUExEREREREWnESjURERERERGRRqxUExEREREREWnESjURERERERGRRqxUExEREREREWn0/wFmXVzqh5YH6wAAAABJRU5ErkJggg==",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.figure()\n",
+ "\n",
+ "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
+ "\n",
+ "data_in.plot(ax = axes[0])\n",
+ "axes[0].set_title(\"ds_in\")\n",
+ "\n",
+ "data_out.plot(ax = axes[1])\n",
+ "axes[1].set_title(\"ds_out regridded\")\n",
+ "\n",
+ "plt.tight_layout();"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "585f9690-e5fb-4383-97f4-a0d37ac2fddb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def compute_integral(data, ds, R_earth):\n",
+ "\n",
+ " lon_b = np.deg2rad(ds['lon_b'].values)\n",
+ " lat_b = np.deg2rad(ds['lat_b'].values)\n",
+ " area = np.zeros((ds.sizes['y'], ds.sizes['x']))\n",
+ " for j in range(ds.sizes['y']):\n",
+ " for i in range(ds.sizes['x']):\n",
+ " dlon = lon_b[j, i+1] - lon_b[j, i]\n",
+ " dlat = np.sin(lat_b[j+1, i]) - np.sin(lat_b[j, i])\n",
+ " area[j, i] = R_earth**2 * dlon * dlat\n",
+ "\n",
+ " return np.sum(data.values * area)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "6424fb9b-b096-427e-86e2-d828bbec001d",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "0.0\n",
+ "0.0\n",
+ "510064471909788.25\n",
+ "510064471909788.25\n"
+ ]
+ }
+ ],
+ "source": [
+ "R_earth = 6.371e6\n",
+ "\n",
+ "# test that compute_integral works;\n",
+ "# gives ones as input and check if it gives the area of Earth\n",
+ "\n",
+ "print(compute_integral(xr.ones_like(data_in, dtype=float), ds_in, R_earth) - 4*np.pi*R_earth**2)\n",
+ "print(compute_integral(xr.ones_like(data_out, dtype=float), ds_out, R_earth) - 4*np.pi*R_earth**2)\n",
+ "print(compute_integral(xr.ones_like(data_in, dtype=float), ds_in, R_earth))\n",
+ "print(compute_integral(xr.ones_like(data_out, dtype=float), ds_out, R_earth))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "48d62d37-4ee1-4bfd-9a56-879dec07b09b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "integral_in = compute_integral(data_in, ds_in, R_earth)\n",
+ "integral_out = compute_integral(data_out, ds_out, R_earth)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "81ab1c27-ffca-47c7-854b-ff81445e92c5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "integral ds_in : -11089929769375.488\n",
+ "integral ds_out: -11090665232069.717\n",
+ "abs error: 735462694.2285156\n",
+ "rel error: 6.631806598626747e-05\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"integral ds_in : \", integral_in)\n",
+ "print(\"integral ds_out: \", integral_out)\n",
+ "\n",
+ "print(\"abs error: \", np.abs(integral_in - integral_out))\n",
+ "print(\"rel error: \", np.abs(integral_in - integral_out) / np.abs(integral_in))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1ce6ac14-f9f8-4d3f-9eab-480db7114df5",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.12.8"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/validation/regridding/xesmf_regridding.jl b/validation/regridding/xesmf_regridding.jl
new file mode 100644
index 0000000000..b281341d65
--- /dev/null
+++ b/validation/regridding/xesmf_regridding.jl
@@ -0,0 +1,32 @@
+using Oceananigans
+using XESMF
+
+arch = CPU()
+
+z = (-1, 0)
+
+radius = Oceananigans.defaults.planet_radius
+
+llg_coarse = LatitudeLongitudeGrid(arch; z, radius,
+ size = (180, 90, 1),
+ longitude = (0, 360),
+ latitude = (-90, 90))
+
+llg_fine = LatitudeLongitudeGrid(arch; z, radius,
+ size = (360, 180, 1),
+ longitude = (0, 360),
+ latitude = (-90, 90))
+
+src_field = CenterField(llg_coarse)
+dst_field = CenterField(llg_fine)
+
+λ₀, φ₀ = 150, 30 # degrees
+width = 12 # degrees
+set!(src_field, (λ, φ, z) -> exp(-((λ - λ₀)^2 + (φ - φ₀)^2) / 2width^2) - 2exp(-((λ - 270)^2 + (φ + 20)^2) / 2width^2))
+
+regridder = XESMF.Regridder(dst_field, src_field, method="conservative")
+
+regrid!(dst_field, regridder, src_field)
+
+println("integral src_field = ", first(Field(Integral(src_field, dims=(1, 2)))))
+println("integral dst_field = ", first(Field(Integral(dst_field, dims=(1, 2)))))