chore(release): 6.0.0

# [6.0.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.13.0...v6.0.0) (2026-01-19) ### Bug Fixes * Update Arome to translate latlon to lambert ([655abeb](655abebe52))
Merge branch 'simkir/arome' into 'main'
2026-01-19 14:13:19 +00:00 · 2026-01-19 15:10:29 +01:00 · 2026-01-07 09:34:08 +01:00 · 2026-01-06 10:23:49 +01:00 · 2025-12-17 17:40:12 +01:00 · 2025-11-05 11:09:46 +00:00
47 changed files with 2980 additions and 4120 deletions
--- a/.build/Build.fs
+++ b/.build/Build.fs
@@ -11,23 +11,23 @@ let srcPath  = Path.getFullName "src"
 let testPath = Path.getFullName "test"
 let libPath  = Some srcPath

-let deployPath  = Path.getFullName "deploy"
+let distPath  = Path.getFullName "dist"
 let packPath    = Path.getFullName "packages"
 let versionFile = Path.getFullName ".version"

-Target.create "Clean" (fun _ -> Shell.cleanDir deployPath)
+Target.create "Clean" (fun _ -> Shell.cleanDir distPath)

 Target.create "InstallClient" (fun _ ->
-    run npm "install" "."
+    run bun "install" "."
    run dotnet "tool restore" "."
 )

 Target.create "Bundle" (fun _ ->
-    run dotnet $"publish -c Release -o \"{deployPath}\"" srcPath
+    run dotnet $"publish -c Release -o \"{distPath}\"" srcPath
 )

 Target.create "BundleDebug" (fun _ ->
-    run dotnet $"publish -c Debug -o \"{deployPath}\"" srcPath
+    run dotnet $"publish -c Debug -o \"{distPath}\"" srcPath
 )

 Target.create "Pack" (fun _ ->
--- a/.build/Helpers.fs
+++ b/.build/Helpers.fs
@@ -72,16 +72,16 @@ let createProcess exe arg dir =
    |> CreateProcess.ensureExitCode

 let dotnet = createProcess "dotnet"
-let npm =
-    let npmPath =
-        match ProcessUtils.tryFindFileOnPath "npm" with
+let bun =
+    let bunPath =
+        match ProcessUtils.tryFindFileOnPath "bun" with
        | Some path -> path
        | None ->
-            "npm was not found in path. Please install it and make sure it's available from your path. " +
+            "bun was not found in path. Please install it and make sure it's available from your path. " +
            "See https://safe-stack.github.io/docs/quickstart/#install-pre-requisites for more info"
            |> failwith

-    createProcess npmPath
+    createProcess bunPath

 let run proc arg dir =
    proc arg dir
--- a/.config/dotnet-tools.json
+++ b/.config/dotnet-tools.json
@@ -1,18 +0,0 @@
-{
-  "version": 1,
-  "isRoot": true,
-  "tools": {
-    "fable": {
-      "version": "3.7.0",
-      "commands": [
-        "fable"
-      ]
-    },
-    "fantomas-tool": {
-      "version": "4.6.4",
-      "commands": [
-        "fantomas"
-      ]
-    }
-  }
-}
--- a/.devcontainer/Dockerfile
+++ b/.devcontainer/Dockerfile
@@ -1,27 +1,30 @@
-FROM mcr.microsoft.com/dotnet/sdk:6.0
+FROM mcr.microsoft.com/dotnet/sdk:9.0

-# Add keys and sources lists
-RUN curl -sL https://deb.nodesource.com/setup_14.x | bash
-RUN curl -sS https://dl.yarnpkg.com/debian/pubkey.gpg | apt-key add -
-RUN echo "deb https://dl.yarnpkg.com/debian/ stable main" \
-    | tee /etc/apt/sources.list.d/yarn.list
+# Bun version
+ARG BUN_INSTALL=/usr/local
+ARG BUN_VERSION=bun-v1.2.16

 # Install node, 7zip, yarn, git, process tools
 RUN apt-get update \
-    && apt-get install -y nodejs p7zip-full git procps ssh-client
+    && apt-get install -y p7zip-full git procps ssh-client unzip
+
+# Install Bun
+RUN set -eux; \
+    curl -fsSL https://bun.sh/install > /usr/local/bin/install-bun \
+    && chmod +x /usr/local/bin/install-bun \
+    && /usr/local/bin/install-bun $BUN_VERSION debug-info
+
+ENV BUN_INSTALL=/usr/local

 # Clean up
 RUN apt-get autoremove -y \
    && apt-get clean -y \
    && rm -rf /var/lib/apt/lists/*

-# Install dotnet tools
-RUN dotnet tool install fable -g
-
 # Trouble brewing
 RUN rm /etc/ssl/openssl.cnf

-# add dotnet tools to path to pick up fake and paket installation
+# Add dotnet tools to path to pick up fake and paket installation
 ENV PATH="/root/.dotnet/tools:${PATH}"

 # Copy endpoint specific user settings into container to specify
--- a/.editorconfig
+++ b/.editorconfig
@@ -7,26 +7,30 @@ charset = utf-8
 trim_trailing_whitespace = true
 insert_final_newline = false

+[*.js]
+indent_size = 2
+max_line_length= 80
+
+[*.nix]
+indent_size = 2
+max_line_length= 80
+
 [*.fs]
-max_line_length=120
-# Feliz style
-fsharp_single_argument_web_mode=true
-fsharp_space_before_colon=false
-fsharp_max_if_then_else_short_width=60
-fsharp_max_infix_operator_expression=50
-fsharp_max_record_width=70
-fsharp_max_record_number_of_items=1
-fsharp_max_array_or_list_width=70
-fsharp_max_array_or_list_number_of_items=1
-fsharp_max_value_binding_width=70
-fsharp_max_function_binding_width=40
-fsharp_max_dot_get_expression_width=50
-fsharp_multiline_block_brackets_on_same_column=true
-fsharp_newline_between_type_definition_and_members=false
-fsharp_max_elmish_width=40
-fsharp_align_function_signature_to_indentation=false
-fsharp_alternative_long_member_definitions=false
-fsharp_multi_line_lambda_closing_newline=false
-fsharp_disable_elmish_syntax=false
-fsharp_keep_indent_in_branch=false
-fsharp_blank_lines_around_nested_multiline_expressions=false
+max_line_length= 120
+
+fsharp_max_if_then_else_short_width = 60
+fsharp_max_infix_operator_expression = 80
+
+fsharp_space_before_uppercase_invocation = true
+fsharp_blank_lines_around_nested_multiline_expressions = false
+fsharp_newline_between_type_definition_and_members = false
+fsharp_multiline_bracket_style = stroustrup
+fsharp_multi_line_lambda_closing_newline = true
+
+fsharp_array_or_list_multiline_formatter = character_width
+fsharp_max_array_or_list_width = 70
+fsharp_max_array_or_list_number_of_items = 3
+
+fsharp_record_multiline_formatter = number_of_items
+fsharp_max_record_number_of_items = 3
+fsharp_max_record_width = 70
--- a/.envrc
+++ b/.envrc
@@ -0,0 +1,3 @@
+export NPINS_DIRECTORY="nix"
+
+use_nix
--- a/.gitignore
+++ b/.gitignore
@@ -15,3 +15,6 @@ deploy
 .ionide/
 *.db
 build.fsx.lock
+dist/
+.direnv/
+result*
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,10 +1,11 @@
 variables:
-  DEPLOY_NAME: default
-  DEPLOY_NAMESPACE: default
+  SDK_VERSION: 9.0
  SKIP_TESTS: "true"

 include:
  - project: oceanbox/gitlab-ci
-    ref: main
+    ref: v4.1
    file: DotnetPackage.gitlab-ci.yml
-
+    inputs:
+        project-name: oceanbox.fvcomkit
+        project-dir: .
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,8 +0,0 @@
-[submodule "submodules/ProjNet.FSharp"]
-	path = submodules/ProjNet.FSharp
-	url = ../ProjNet.FSharp
-    shallow = true
-[submodule "submodules/FsKdTree"]
-	path = submodules/FsKdTree
-	url = https://gitlab.com/serit/libs/FsKDTree.git
-    shallow = true
--- a/.version
+++ b/.version
@@ -1 +1 @@
-3.12.2
+6.0.0
--- a/Build.fsproj
+++ b/Build.fsproj
@@ -1,17 +1,17 @@
 <Project Sdk="Microsoft.NET.Sdk">
  <PropertyGroup>
    <OutputType>Exe</OutputType>
-    <TargetFramework>net6.0</TargetFramework>
+    <TargetFramework>net9.0</TargetFramework>
  </PropertyGroup>
  <ItemGroup>
    <Compile Include=".build/Helpers.fs" />
    <Compile Include=".build/Build.fs" />
  </ItemGroup>
  <ItemGroup>
-    <PackageReference Include="Fake.Core.Target" Version="5.23.0" />
-    <PackageReference Include="Fake.DotNet.Cli" Version="5.23.0" />
-    <PackageReference Include="Fake.IO.FileSystem" Version="5.23.0" />
-    <PackageReference Include="Farmer" Version="1.7.7" />
-    <PackageReference Update="FSharp.Core" Version="6.0.4" />
+    <PackageReference Include="Fake.Core.Target" Version="6.1.3" />
+    <PackageReference Include="Fake.DotNet.Cli" Version="6.1.3" />
+    <PackageReference Include="Fake.IO.FileSystem" Version="6.1.3" />
+    <PackageReference Include="Farmer" Version="1.9.11" />
+    <PackageReference Update="FSharp.Core" Version="9.0.201" />
  </ItemGroup>
-</Project>
+</Project>
--- a/Oceanbox.FvcomKit.sln
+++ b/Oceanbox.FvcomKit.sln
@@ -1,77 +0,0 @@
-
-Microsoft Visual Studio Solution File, Format Version 12.00
-# Visual Studio 15
-VisualStudioVersion = 15.0.27004.2005
-MinimumVisualStudioVersion = 15.0.26124.0
-Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{C29C6F32-3A30-4071-9B4A-8FBCAAA5993A}"
-	ProjectSection(SolutionItems) = preProject
-		README.md = README.md
-		LICENSE = LICENSE
-		Dockerfile = Dockerfile
-	EndProjectSection
-EndProject
-Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "Build", "Build.fsproj", "{C6824583-FB68-4F69-8117-6B29637A3B96}"
-EndProject
-Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "src", "src/Oceanbox.FvcomKit.fsproj", "{662A0CDC-7E82-4157-AD25-469DD7ABAA69}"
-EndProject
-Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "test", "test\Tests.fsproj", "{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}"
-EndProject
-Global
-	GlobalSection(SolutionConfigurationPlatforms) = preSolution
-		Debug|Any CPU = Debug|Any CPU
-		Debug|x64 = Debug|x64
-		Debug|x86 = Debug|x86
-		Release|Any CPU = Release|Any CPU
-		Release|x64 = Release|x64
-		Release|x86 = Release|x86
-	EndGlobalSection
-	GlobalSection(SolutionProperties) = preSolution
-		HideSolutionNode = FALSE
-	EndGlobalSection
-	GlobalSection(ProjectConfigurationPlatforms) = postSolution
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Debug|Any CPU.Build.0 = Debug|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Debug|x64.ActiveCfg = Debug|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Debug|x64.Build.0 = Debug|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Debug|x86.ActiveCfg = Debug|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Debug|x86.Build.0 = Debug|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Release|Any CPU.ActiveCfg = Release|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Release|Any CPU.Build.0 = Release|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Release|x64.ActiveCfg = Release|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Release|x64.Build.0 = Release|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Release|x86.ActiveCfg = Release|Any CPU
-		{C6824583-FB68-4F69-8117-6B29637A3B96}.Release|x86.Build.0 = Release|Any CPU
-
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Debug|Any CPU.Build.0 = Debug|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Debug|x64.ActiveCfg = Debug|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Debug|x64.Build.0 = Debug|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Debug|x86.ActiveCfg = Debug|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Debug|x86.Build.0 = Debug|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Release|Any CPU.ActiveCfg = Release|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Release|Any CPU.Build.0 = Release|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Release|x64.ActiveCfg = Release|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Release|x64.Build.0 = Release|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Release|x86.ActiveCfg = Release|Any CPU
-		{662A0CDC-7E82-4157-AD25-469DD7ABAA69}.Release|x86.Build.0 = Release|Any CPU
-
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Debug|Any CPU.Build.0 = Debug|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Debug|x64.ActiveCfg = Debug|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Debug|x64.Build.0 = Debug|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Debug|x86.ActiveCfg = Debug|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Debug|x86.Build.0 = Debug|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Release|Any CPU.ActiveCfg = Release|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Release|Any CPU.Build.0 = Release|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Release|x64.ActiveCfg = Release|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Release|x64.Build.0 = Release|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Release|x86.ActiveCfg = Release|Any CPU
-		{BCBD73E2-7170-4A85-BFD8-B76F056D4D49}.Release|x86.Build.0 = Release|Any CPU
-
-	EndGlobalSection
-	GlobalSection(NestedProjects) = preSolution
-	EndGlobalSection
-	GlobalSection(ExtensibilityGlobals) = postSolution
-		SolutionGuid = {79A6998D-BCE6-4EC5-ADBC-69234C0D2EC5}
-	EndGlobalSection
-EndGlobal
--- a/Oceanbox.FvcomKit.slnx
+++ b/Oceanbox.FvcomKit.slnx
@@ -0,0 +1,16 @@
+<Solution>
+  <Configurations>
+    <Platform Name="Any CPU" />
+    <Platform Name="x64" />
+    <Platform Name="x86" />
+  </Configurations>
+  <Folder Name="/Solution Items/">
+    <File Path="LICENSE" />
+    <File Path="README.md" />
+    <File Path="shell.nix" />
+  </Folder>
+  <Folder Name="/submodules/" />
+  <Project Path="Build.fsproj" />
+  <Project Path="src/Oceanbox.FvcomKit.fsproj" />
+  <Project Path="xtest/xtest.fsproj" />
+</Solution>
--- a/RELEASE_NOTES.md
+++ b/RELEASE_NOTES.md
@@ -1,5 +1,335 @@
 # Changelog

+# [6.0.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.13.0...v6.0.0) (2026-01-19)
+
+
+### Bug Fixes
+
+* Update Arome to translate latlon to lambert ([655abeb](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/655abebe526a3587c26429e41130d78f326fa524))
+
+# [5.13.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.12.2...v5.13.0) (2025-11-05)
+
+
+### Features
+
+* Add ElemsAroundElem to NeighborIndex ([71e861e](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/71e861e4174a57d11cfd85ac04f69c51fc985b4c))
+
+## [5.12.2](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.12.1...v5.12.2) (2025-09-01)
+
+
+### Bug Fixes
+
+* Bump Oceanbox.SDSLite to 2.8.0 and use bun for SR ([0c55b2d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/0c55b2dcbdb89337abdbe430a30461d09dbb0bd2))
+
+## [5.12.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.12.0...v5.12.1) (2025-05-02)
+
+
+### Bug Fixes
+
+* include edge point in isInsideTriangle ([c89f35b](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/c89f35bc6e316ac19d2de4ef35dbb78302881373))
+* remove unused double precision function ([8715f4d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/8715f4d8c32d04f6d1c77502435234dadd0863ae))
+
+# [5.12.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.11.0...v5.12.0) (2025-03-07)
+
+
+### Features
+
+* add Grid.toLonLat function ([06a4aea](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/06a4aeabf1bcd2824db74a1851f0703797fdacf8))
+
+# [5.11.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.10.0...v5.11.0) (2025-03-06)
+
+
+### Bug Fixes
+
+* init sha once explicitly from sha1 byte[] ([33b7b99](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/33b7b999c8aaff3190dbe8a5e4d3f2e8e1cbf15e))
+
+
+### Features
+
+* rename grid sha1 to hash for better generality ([bfbaa3a](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/bfbaa3aeff66b39526cf817009729df3fe4966cc))
+
+# [5.10.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.9.1...v5.10.0) (2025-03-06)
+
+
+### Features
+
+* init sha once explicitly from bingrid ([2e87294](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/2e87294c46579b35096c80b50bc16b89a1352b93))
+
+## [5.9.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.9.0...v5.9.1) (2025-03-06)
+
+
+### Bug Fixes
+
+* use sha from bingrid if it exists ([abdb949](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/abdb94911343c7dff214047ec8bacf740aeffe97))
+
+# [5.9.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.8.0...v5.9.0) (2025-03-06)
+
+
+### Features
+
+* add ToGrid() method to extended grid ([21e84d0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/21e84d07cd4cf1593ae3e9c9f65718d221adcf1e))
+
+# [5.8.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.7.0...v5.8.0) (2025-03-06)
+
+
+### Features
+
+* compute sha1 checksum of extended grids ([bf7ae58](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/bf7ae5889b48efb871c8c3b5ee4f07ae6ce34865))
+
+# [5.7.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.6.0...v5.7.0) (2025-02-12)
+
+
+### Features
+
+* add read omega block ([9abd9d4](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/9abd9d45bc3ce5ec2708f9e2dcbb1366190aad9c))
+
+# [5.6.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.5.5...v5.6.0) (2024-11-27)
+
+
+### Features
+
+* update to net9.0 ([17685ac](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/17685ac2cd781a05da5064d344ab6e75597d48f9))
+
+## [5.5.5](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.5.4...v5.5.5) (2024-11-27)
+
+
+### Bug Fixes
+
+* add velocity to nodal function ([ed0ac79](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/ed0ac797d378c7acc0dba29cc5d95f10d1e0fe3a))
+
+## [5.5.4](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.5.3...v5.5.4) (2024-05-01)
+
+
+### Bug Fixes
+
+* Norshelf filenames changed due to changes in Thredds ([0f5de91](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/0f5de91d39e7fae82b17bfc5fd014f731c7d2b9d))
+
+## [5.5.3](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.5.2...v5.5.3) (2024-02-23)
+
+
+### Bug Fixes
+
+* bug in isInsideTriangle ([8d7ab41](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/8d7ab4164c2817906c539b1c632c9ea820de7e5d))
+
+## [5.5.2](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.5.1...v5.5.2) (2024-01-04)
+
+
+### Bug Fixes
+
+* fix ci/cd deploy name ([35e1f43](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/35e1f43fc00dd31c2ec442a3e81d6b715011e66e))
+* update ci/cd setup to v2 ([3f61f1d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/3f61f1d1a9c98dfc6be98540cfbc65e90b50cc36))
+
+## [5.5.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.5.0...v5.5.1) (2024-01-03)
+
+
+### Bug Fixes
+
+* bump to sdk_8 in shell.nix ([bc13065](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/bc130658e97e7811c26aaa132d945c4a0258fa9c))
+* update package sdslite.oceanbox-2.7.3 ([aba0917](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/aba0917557b32ad29aabba0d9d59dd8925477a8c))
+
+# [5.5.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.4.1...v5.5.0) (2023-12-31)
+
+
+### Bug Fixes
+
+* update devcontainer to net8.0 ([d75db35](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/d75db35d8bba9283091447ae08ccf9ec07213f7a))
+
+
+### Features
+
+* upgrade to net8.0 ([cfbf2b1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/cfbf2b1d61e1a68c4516ffe36c0ac28526fb88e1))
+
+## [5.4.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.4.0...v5.4.1) (2023-10-04)
+
+
+### Bug Fixes
+
+* fix arome variable names ([0f2b696](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/0f2b69609c5ac8d5cc1e8ddf5e2cdb8e0e2b7ef1))
+
+# [5.4.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.3.0...v5.4.0) (2023-09-25)
+
+
+### Bug Fixes
+
+* bug in bbox center ([90b50c3](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/90b50c3c8359536ffffb8737cde30682f079fbd3))
+* remove depricated grid.projection ([f38f06c](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/f38f06c936d92ca2598f16b56e3d7025d8efb4ae))
+* remove print ([d2e362d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/d2e362df2019f70d75a3d9a79b7f0fee3992e574))
+* rescale/translate grid ([5f5ad1b](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/5f5ad1b9e41fbb502db41b9fe327ec970a1a608e))
+* update ProjNet (error on unknown proj) ([5479c8b](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/5479c8b598c2c1474f86d841fad1fac8e6b431d5))
+
+
+### Features
+
+* read tauc from fvcom file ([4ca7fff](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/4ca7fff09ea4eea3a8af66bd71907de96df642c3))
+
+# [5.3.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.2.0...v5.3.0) (2023-09-06)
+
+
+### Features
+
+* add functtions to read grid in lon-lat format ([902ac08](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/902ac080e9b74167825fbccb086f99f6bad203c2))
+
+# [5.2.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.1.1...v5.2.0) (2023-08-28)
+
+
+### Bug Fixes
+
+* getTimeInDays returns single, not float ([434ab6d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/434ab6d5f43390e75f9760bd9310bbbeed678975))
+* remove FsKdTree references for faster builds ([86384bf](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/86384bfd9f64e768c93cda2b04c2787897d2b20f))
+* update deps ([460f8f1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/460f8f10eb3aed704c59eaddf587268b43373896))
+* update deps ([bb57404](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/bb5740452375984cd97139d64a8a3e5c204f6661))
+
+
+### Features
+
+* add getTimeInDays to Fvcom accessors ([302c4f1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/302c4f10488cf017501427e044b6f3669fc28258))
+* add getTimsSpanSinceStart function ([764bded](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/764bdedb453c4a8e2bbec3c1cd43ed314fa08e94))
+* add Singlular module for time-series (later add Plural) ([eb2a981](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/eb2a98157653d8a7faadbcc7c6106de8befd867a))
+* read data in sigma layer blocks ([e7bb554](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/e7bb5540ef74a88efe259ef8063c46ec33af7420))
+* read necessary grid prop from file ([ff947bb](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/ff947bb3644af7d09039df5898ce825ad1ca88bf))
+* remove projection from grid(s), now in archmesiter ([36bd52d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/36bd52d3f174ea76f155c680066dca48e9d18a12))
+* upgrade to net7.0 ([2be8397](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/2be8397297e586ca7292c40dc65e0484a44483d3))
+* working arome querying ([4556a9f](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/4556a9fd9a7f8bd9f15719577336a894e2435ba8))
+
+## [5.1.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.1.0...v5.1.1) (2023-02-01)
+
+
+### Bug Fixes
+
+* compute and not read center variables ([a6a8828](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/a6a8828ae9f403469d827f767cd001f7191370b5))
+* compute and not read siglay_center ([a9d7662](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/a9d7662a2b131bb43740d8a279e4162db5a99894))
+
+# [5.1.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.0.3...v5.1.0) (2023-01-16)
+
+
+### Bug Fixes
+
+* update node in devcontainer ([347a730](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/347a730e135f407d2635a29f63a25796ebf6b2fb))
+
+
+### Features
+
+* add readOmega for vertical sigma velocity ([7d81ac7](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/7d81ac771a0b14ba5a3da97e4986610056969ff9))
+
+## [5.0.3](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.0.2...v5.0.3) (2023-01-02)
+
+
+### Bug Fixes
+
+* using S_rho in readVerticalGrid ([886234b](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/886234b1659a214755b8e439e99d003033f992fd))
+
+## [5.0.2](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.0.1...v5.0.2) (2023-01-02)
+
+
+### Bug Fixes
+
+* use S_rho as vertical roms coordinate ([f851e54](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/f851e5411b5c248385be216de1409fda68603a7b))
+
+## [5.0.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v5.0.0...v5.0.1) (2022-12-19)
+
+
+### Bug Fixes
+
+* update ProjNet.FSharp ([0bf255d](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/0bf255d51056bfe506d3c9da7f2765d74aa75244))
+
+# [5.0.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.6.0...v5.0.0) (2022-12-07)
+
+
+### Bug Fixes
+
+* update ProjNet.FSharp ([ab63c39](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/ab63c39993906dbdf80f9fd19ec636981a348d0f))
+
+# [4.6.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.5.3...v4.6.0) (2022-12-02)
+
+
+### Features
+
+* add tryGetNode grid method ([8c1c3e7](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/8c1c3e750f67d3a4c50b0a3bf1f74a8bfa792410))
+
+## [4.5.3](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.5.2...v4.5.3) (2022-12-02)
+
+
+### Bug Fixes
+
+* make siglev/lay readers faster ([a2b6a3b](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/a2b6a3b04fa71b38beffa7c06fcdc5ae076b3299))
+
+## [4.5.2](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.5.1...v4.5.2) (2022-11-16)
+
+
+### Bug Fixes
+
+* fix obc indexing off by one ([0bf1889](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/0bf18891b26631c2a1065cf49415e59a96efe4fa))
+
+## [4.5.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.5.0...v4.5.1) (2022-11-15)
+
+
+### Bug Fixes
+
+* update ProjNet.FSharp ([1939213](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/1939213ec40f0e3ca46676e97d3e9a05a42d92d1))
+
+# [4.5.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.4.1...v4.5.0) (2022-11-03)
+
+
+### Features
+
+* add functions to ger number of sigmas ([7c977ab](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/7c977abaa3b7a36fade925db3aaa99122ffdb482))
+
+## [4.4.1](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.4.0...v4.4.1) (2022-10-18)
+
+
+### Bug Fixes
+
+* store Fvcom field as single, but use as float ([8cd202a](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/8cd202a6993a7d59a138d4ccefe8251bf47dc0b2))
+
+# [4.4.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.3.0...v4.4.0) (2022-10-17)
+
+
+### Features
+
+* add readUV(W) ranges ([e2c977b](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/e2c977b843552f573cb672b60c2aa78ae2db8732))
+
+# [4.3.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.2.0...v4.3.0) (2022-10-12)
+
+
+### Features
+
+* compute circumscribed circle around element ([3c9a7a7](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/3c9a7a701ef1cdc3d168608f4d85dc2b08a8de6e))
+
+# [4.2.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.1.0...v4.2.0) (2022-10-05)
+
+
+### Features
+
+* add Fvcom.getTime ([cbe31b9](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/cbe31b94b50a6cf7a066e095af27ff073855c119))
+
+# [4.1.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v4.0.0...v4.1.0) (2022-09-30)
+
+
+### Features
+
+* project on boundary ([0a8b31c](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/0a8b31c1f3b85ccdf9b5e108499352aabe7e0e79))
+
+# [4.0.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v3.14.0...v4.0.0) (2022-09-17)
+
+
+### Bug Fixes
+
+* reformat src with fantoms ([7955a76](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/7955a76178cbe890901cf825bf9fc06013cf56ee))
+
+# [3.14.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v3.13.0...v3.14.0) (2022-09-12)
+
+
+### Features
+
+* add readUVW(s) and fix bug in readUVs ([80b2651](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/80b2651f2441b8b85037995a400c38db581c18bd))
+
+# [3.13.0](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v3.12.2...v3.13.0) (2022-09-09)
+
+
+### Features
+
+* add elemental accessors for U/V and WW ([060c775](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/commit/060c7753a2f04cf6388a9bec62e92fe9759de6e2))
+
 ## [3.12.2](https://gitlab.com/oceanbox/Oceanbox.FvcomKit/compare/v3.12.1...v3.12.2) (2022-09-09)


--- a/bun.lock
+++ b/bun.lock
@@ -0,0 +1,302 @@
+{
+  "lockfileVersion": 1,
+  "configVersion": 0,
+  "workspaces": {
+    "": {
+      "devDependencies": {
+        "@semantic-release/changelog": "^6.0.1",
+        "@semantic-release/exec": "^6.0.3",
+        "@semantic-release/git": "^10.0.1",
+        "@semantic-release/gitlab": "^7.0.4",
+        "semantic-release-dotnet": "^1.0.0",
+      },
+    },
+  },
+  "packages": {
+    "@babel/code-frame": ["@babel/code-frame@7.16.7", "", { "dependencies": { "@babel/highlight": "^7.16.7" } }, "sha512-iAXqUn8IIeBTNd72xsFlgaXHkMBMt6y4HJp1tIaK465CWLT/fG1aqB7ykr95gHHmlBdGbFeWWfyB4NJJ0nmeIg=="],
+
+    "@babel/helper-validator-identifier": ["@babel/helper-validator-identifier@7.16.7", "", {}, "sha512-hsEnFemeiW4D08A5gUAZxLBTXpZ39P+a+DGDsHw1yxqyQ/jzFEnxf5uTEGp+3bzAbNOxU1paTgYS4ECU/IgfDw=="],
+
+    "@babel/highlight": ["@babel/highlight@7.16.7", "", { "dependencies": { "@babel/helper-validator-identifier": "^7.16.7", "chalk": "^2.0.0", "js-tokens": "^4.0.0" } }, "sha512-aKpPMfLvGO3Q97V0qhw/V2SWNWlwfJknuwAunU7wZLSfrM4xTBvg7E5opUVi1kJTBKihE38CPg4nBiqX83PWYw=="],
+
+    "@nodelib/fs.scandir": ["@nodelib/fs.scandir@2.1.5", "", { "dependencies": { "@nodelib/fs.stat": "2.0.5", "run-parallel": "^1.1.9" } }, "sha512-vq24Bq3ym5HEQm2NKCr3yXDwjc7vTsEThRDnkp2DK9p1uqLR+DHurm/NOTo0KG7HYHU7eppKZj3MyqYuMBf62g=="],
+
+    "@nodelib/fs.stat": ["@nodelib/fs.stat@2.0.5", "", {}, "sha512-RkhPPp2zrqDAQA/2jNhnztcPAlv64XdhIp7a7454A5ovI7Bukxgt7MX7udwAu3zg1DcpPU0rz3VV1SeaqvY4+A=="],
+
+    "@nodelib/fs.walk": ["@nodelib/fs.walk@1.2.8", "", { "dependencies": { "@nodelib/fs.scandir": "2.1.5", "fastq": "^1.6.0" } }, "sha512-oGB+UxlgWcgQkgwo8GcEGwemoTFt3FIO9ababBmaGwXIoBKZ+GTy0pP185beGg7Llih/NSHSV2XAs1lnznocSg=="],
+
+    "@semantic-release/changelog": ["@semantic-release/changelog@6.0.1", "", { "dependencies": { "@semantic-release/error": "^3.0.0", "aggregate-error": "^3.0.0", "fs-extra": "^9.0.0", "lodash": "^4.17.4" } }, "sha512-FT+tAGdWHr0RCM3EpWegWnvXJ05LQtBkQUaQRIExONoXjVjLuOILNm4DEKNaV+GAQyJjbLRVs57ti//GypH6PA=="],
+
+    "@semantic-release/error": ["@semantic-release/error@3.0.0", "", {}, "sha512-5hiM4Un+tpl4cKw3lV4UgzJj+SmfNIDCLLw0TepzQxz9ZGV5ixnqkzIVF+3tp0ZHgcMKE+VNGHJjEeyFG2dcSw=="],
+
+    "@semantic-release/exec": ["@semantic-release/exec@6.0.3", "", { "dependencies": { "@semantic-release/error": "^3.0.0", "aggregate-error": "^3.0.0", "debug": "^4.0.0", "execa": "^5.0.0", "lodash": "^4.17.4", "parse-json": "^5.0.0" } }, "sha512-bxAq8vLOw76aV89vxxICecEa8jfaWwYITw6X74zzlO0mc/Bgieqx9kBRz9z96pHectiTAtsCwsQcUyLYWnp3VQ=="],
+
+    "@semantic-release/git": ["@semantic-release/git@10.0.1", "", { "dependencies": { "@semantic-release/error": "^3.0.0", "aggregate-error": "^3.0.0", "debug": "^4.0.0", "dir-glob": "^3.0.0", "execa": "^5.0.0", "lodash": "^4.17.4", "micromatch": "^4.0.0", "p-reduce": "^2.0.0" } }, "sha512-eWrx5KguUcU2wUPaO6sfvZI0wPafUKAMNC18aXY4EnNcrZL86dEmpNVnC9uMpGZkmZJ9EfCVJBQx4pV4EMGT1w=="],
+
+    "@semantic-release/gitlab": ["@semantic-release/gitlab@7.0.4", "", { "dependencies": { "@semantic-release/error": "^3.0.0", "aggregate-error": "^3.0.0", "debug": "^4.0.0", "dir-glob": "^3.0.0", "escape-string-regexp": "^3.0.0", "form-data": "^4.0.0", "fs-extra": "^10.0.0", "globby": "^11.0.0", "got": "^11.0.0", "lodash": "^4.17.11", "parse-path": "^4.0.0", "url-join": "^4.0.0" } }, "sha512-TL6kT526+ir/uehMFdTlJNXUj+p+SjPAYUkit6lh5Rs8kxeHQ01bgmpYLQlc94ZDpy9x2Tzcb/NRwKojkmLG4A=="],
+
+    "@sindresorhus/is": ["@sindresorhus/is@4.3.0", "", {}, "sha512-wwOvh0eO3PiTEivGJWiZ+b946SlMSb4pe+y+Ur/4S87cwo09pYi+FWHHnbrM3W9W7cBYKDqQXcrFYjYUCOJUEQ=="],
+
+    "@szmarczak/http-timer": ["@szmarczak/http-timer@4.0.6", "", { "dependencies": { "defer-to-connect": "^2.0.0" } }, "sha512-4BAffykYOgO+5nzBWYwE3W90sBgLJoUPRWWcL8wlyiM8IB8ipJz3UMJ9KXQd1RKQXpKp8Tutn80HZtWsu2u76w=="],
+
+    "@types/cacheable-request": ["@types/cacheable-request@6.0.2", "", { "dependencies": { "@types/http-cache-semantics": "*", "@types/keyv": "*", "@types/node": "*", "@types/responselike": "*" } }, "sha512-B3xVo+dlKM6nnKTcmm5ZtY/OL8bOAOd2Olee9M1zft65ox50OzjEHW91sDiU9j6cvW8Ejg1/Qkf4xd2kugApUA=="],
+
+    "@types/glob": ["@types/glob@7.2.0", "", { "dependencies": { "@types/minimatch": "*", "@types/node": "*" } }, "sha512-ZUxbzKl0IfJILTS6t7ip5fQQM/J3TJYubDm3nMbgubNNYS62eXeUpoLUC8/7fJNiFYHTrGPQn7hspDUzIHX3UA=="],
+
+    "@types/http-cache-semantics": ["@types/http-cache-semantics@4.0.1", "", {}, "sha512-SZs7ekbP8CN0txVG2xVRH6EgKmEm31BOxA07vkFaETzZz1xh+cbt8BcI0slpymvwhx5dlFnQG2rTlPVQn+iRPQ=="],
+
+    "@types/keyv": ["@types/keyv@3.1.3", "", { "dependencies": { "@types/node": "*" } }, "sha512-FXCJgyyN3ivVgRoml4h94G/p3kY+u/B86La+QptcqJaWtBWtmc6TtkNfS40n9bIvyLteHh7zXOtgbobORKPbDg=="],
+
+    "@types/minimatch": ["@types/minimatch@3.0.5", "", {}, "sha512-Klz949h02Gz2uZCMGwDUSDS1YBlTdDDgbWHi+81l29tQALUtvz4rAYi5uoVhE5Lagoq6DeqAUlbrHvW/mXDgdQ=="],
+
+    "@types/node": ["@types/node@17.0.9", "", {}, "sha512-5dNBXu/FOER+EXnyah7rn8xlNrfMOQb/qXnw4NQgLkCygKBKhdmF/CA5oXVOKZLBEahw8s2WP9LxIcN/oDDRgQ=="],
+
+    "@types/responselike": ["@types/responselike@1.0.0", "", { "dependencies": { "@types/node": "*" } }, "sha512-85Y2BjiufFzaMIlvJDvTTB8Fxl2xfLo4HgmHzVBz08w4wDePCTjYw66PdrolO0kzli3yam/YCgRufyo1DdQVTA=="],
+
+    "aggregate-error": ["aggregate-error@3.1.0", "", { "dependencies": { "clean-stack": "^2.0.0", "indent-string": "^4.0.0" } }, "sha512-4I7Td01quW/RpocfNayFdFVk1qSuoh0E7JrbRJ16nH01HhKFQ88INq9Sd+nd72zqRySlr9BmDA8xlEJ6vJMrYA=="],
+
+    "ansi-styles": ["ansi-styles@3.2.1", "", { "dependencies": { "color-convert": "^1.9.0" } }, "sha512-VT0ZI6kZRdTh8YyJw3SMbYm/u+NqfsAxEpWO0Pf9sq8/e94WxxOpPKx9FR1FlyCtOVDNOQ+8ntlqFxiRc+r5qA=="],
+
+    "array-union": ["array-union@2.1.0", "", {}, "sha512-HGyxoOTYUyCM6stUe6EJgnd4EoewAI7zMdfqO+kGjnlZmBDz/cR5pf8r/cR4Wq60sL/p0IkcjUEEPwS3GFrIyw=="],
+
+    "asynckit": ["asynckit@0.4.0", "", {}, "sha1-x57Zf380y48robyXkLzDZkdLS3k="],
+
+    "at-least-node": ["at-least-node@1.0.0", "", {}, "sha512-+q/t7Ekv1EDY2l6Gda6LLiX14rU9TV20Wa3ofeQmwPFZbOMo9DXrLbOjFaaclkXKWidIaopwAObQDqwWtGUjqg=="],
+
+    "balanced-match": ["balanced-match@1.0.2", "", {}, "sha512-3oSeUO0TMV67hN1AmbXsK4yaqU7tjiHlbxRDZOpH0KW9+CeX4bRAaX0Anxt0tx2MrpRpWwQaPwIlISEJhYU5Pw=="],
+
+    "brace-expansion": ["brace-expansion@1.1.11", "", { "dependencies": { "balanced-match": "^1.0.0", "concat-map": "0.0.1" } }, "sha512-iCuPHDFgrHX7H2vEI/5xpz07zSHB00TpugqhmYtVmMO6518mCuRMoOYFldEBl0g187ufozdaHgWKcYFb61qGiA=="],
+
+    "braces": ["braces@3.0.2", "", { "dependencies": { "fill-range": "^7.0.1" } }, "sha512-b8um+L1RzM3WDSzvhm6gIz1yfTbBt6YTlcEKAvsmqCZZFw46z626lVj9j1yEPW33H5H+lBQpZMP1k8l+78Ha0A=="],
+
+    "cacheable-lookup": ["cacheable-lookup@5.0.4", "", {}, "sha512-2/kNscPhpcxrOigMZzbiWF7dz8ilhb/nIHU3EyZiXWXpeq/au8qJ8VhdftMkty3n7Gj6HIGalQG8oiBNB3AJgA=="],
+
+    "cacheable-request": ["cacheable-request@7.0.2", "", { "dependencies": { "clone-response": "^1.0.2", "get-stream": "^5.1.0", "http-cache-semantics": "^4.0.0", "keyv": "^4.0.0", "lowercase-keys": "^2.0.0", "normalize-url": "^6.0.1", "responselike": "^2.0.0" } }, "sha512-pouW8/FmiPQbuGpkXQ9BAPv/Mo5xDGANgSNXzTzJ8DrKGuXOssM4wIQRjfanNRh3Yu5cfYPvcorqbhg2KIJtew=="],
+
+    "call-bind": ["call-bind@1.0.2", "", { "dependencies": { "function-bind": "^1.1.1", "get-intrinsic": "^1.0.2" } }, "sha512-7O+FbCihrB5WGbFYesctwmTKae6rOiIzmz1icreWJ+0aA7LJfuqhEso2T9ncpcFtzMQtzXf2QGGueWJGTYsqrA=="],
+
+    "chalk": ["chalk@2.4.2", "", { "dependencies": { "ansi-styles": "^3.2.1", "escape-string-regexp": "^1.0.5", "supports-color": "^5.3.0" } }, "sha512-Mti+f9lpJNcwF4tWV8/OrTTtF1gZi+f8FqlyAdouralcFWFQWF2+NgCHShjkCb+IFBLq9buZwE1xckQU4peSuQ=="],
+
+    "clean-stack": ["clean-stack@2.2.0", "", {}, "sha512-4diC9HaTE+KRAMWhDhrGOECgWZxoevMc5TlkObMqNSsVU62PYzXZ/SMTjzyGAFF1YusgxGcSWTEXBhp0CPwQ1A=="],
+
+    "clone-response": ["clone-response@1.0.2", "", { "dependencies": { "mimic-response": "^1.0.0" } }, "sha1-0dyXOSAxTfZ/vrlCI7TuNQI56Ws="],
+
+    "color-convert": ["color-convert@1.9.3", "", { "dependencies": { "color-name": "1.1.3" } }, "sha512-QfAUtd+vFdAtFQcC8CCyYt1fYWxSqAiK2cSD6zDB8N3cpsEBAvRxp9zOGg6G/SHHJYAT88/az/IuDGALsNVbGg=="],
+
+    "color-name": ["color-name@1.1.3", "", {}, "sha1-p9BVi9icQveV3UIyj3QIMcpTvCU="],
+
+    "combined-stream": ["combined-stream@1.0.8", "", { "dependencies": { "delayed-stream": "~1.0.0" } }, "sha512-FQN4MRfuJeHf7cBbBMJFXhKSDq+2kAArBlmRBvcvFE5BB1HZKXtSFASDhdlz9zOYwxh8lDdnvmMOe/+5cdoEdg=="],
+
+    "concat-map": ["concat-map@0.0.1", "", {}, "sha1-2Klr13/Wjfd5OnMDajug1UBdR3s="],
+
+    "cross-spawn": ["cross-spawn@7.0.3", "", { "dependencies": { "path-key": "^3.1.0", "shebang-command": "^2.0.0", "which": "^2.0.1" } }, "sha512-iRDPJKUPVEND7dHPO8rkbOnPpyDygcDFtWjpeWNCgy8WP2rXcxXL8TskReQl6OrB2G7+UJrags1q15Fudc7G6w=="],
+
+    "debug": ["debug@4.3.3", "", { "dependencies": { "ms": "2.1.2" } }, "sha512-/zxw5+vh1Tfv+4Qn7a5nsbcJKPaSvCDhojn6FEl9vupwK2VCSDtEiEtqr8DFtzYFOdz63LBkxec7DYuc2jon6Q=="],
+
+    "decode-uri-component": ["decode-uri-component@0.2.0", "", {}, "sha1-6zkTMzRYd1y4TNGh+uBiEGu4dUU="],
+
+    "decompress-response": ["decompress-response@6.0.0", "", { "dependencies": { "mimic-response": "^3.1.0" } }, "sha512-aW35yZM6Bb/4oJlZncMH2LCoZtJXTRxES17vE3hoRiowU2kWHaJKFkSBDnDR+cm9J+9QhXmREyIfv0pji9ejCQ=="],
+
+    "defer-to-connect": ["defer-to-connect@2.0.1", "", {}, "sha512-4tvttepXG1VaYGrRibk5EwJd1t4udunSOVMdLSAL6mId1ix438oPwPZMALY41FCijukO1L0twNcGsdzS7dHgDg=="],
+
+    "delayed-stream": ["delayed-stream@1.0.0", "", {}, "sha1-3zrhmayt+31ECqrgsp4icrJOxhk="],
+
+    "dir-glob": ["dir-glob@3.0.1", "", { "dependencies": { "path-type": "^4.0.0" } }, "sha512-WkrWp9GR4KXfKGYzOLmTuGVi1UWFfws377n9cc55/tb6DuqyF6pcQ5AbiHEshaDpY9v6oaSr2XCDidGmMwdzIA=="],
+
+    "end-of-stream": ["end-of-stream@1.4.4", "", { "dependencies": { "once": "^1.4.0" } }, "sha512-+uw1inIHVPQoaVuHzRyXd21icM+cnt4CzD5rW+NC1wjOUSTOs+Te7FOv7AhN7vS9x/oIyhLP5PR1H+phQAHu5Q=="],
+
+    "error-ex": ["error-ex@1.3.2", "", { "dependencies": { "is-arrayish": "^0.2.1" } }, "sha512-7dFHNmqeFSEt2ZBsCriorKnn3Z2pj+fd9kmI6QoWw4//DL+icEBfc0U7qJCisqrTsKTjw4fNFy2pW9OqStD84g=="],
+
+    "escape-string-regexp": ["escape-string-regexp@3.0.0", "", {}, "sha512-11dXIUC3umvzEViLP117d0KN6LJzZxh5+9F4E/7WLAAw7GrHk8NpUR+g9iJi/pe9C0py4F8rs0hreyRCwlAuZg=="],
+
+    "execa": ["execa@5.1.1", "", { "dependencies": { "cross-spawn": "^7.0.3", "get-stream": "^6.0.0", "human-signals": "^2.1.0", "is-stream": "^2.0.0", "merge-stream": "^2.0.0", "npm-run-path": "^4.0.1", "onetime": "^5.1.2", "signal-exit": "^3.0.3", "strip-final-newline": "^2.0.0" } }, "sha512-8uSpZZocAZRBAPIEINJj3Lo9HyGitllczc27Eh5YYojjMFMn8yHMDMaUHE2Jqfq05D/wucwI4JGURyXt1vchyg=="],
+
+    "fast-glob": ["fast-glob@3.2.11", "", { "dependencies": { "@nodelib/fs.stat": "^2.0.2", "@nodelib/fs.walk": "^1.2.3", "glob-parent": "^5.1.2", "merge2": "^1.3.0", "micromatch": "^4.0.4" } }, "sha512-xrO3+1bxSo3ZVHAnqzyuewYT6aMFHRAd4Kcs92MAonjwQZLsK9d0SF1IyQ3k5PoirxTW0Oe/RqFgMQ6TcNE5Ew=="],
+
+    "fastq": ["fastq@1.13.0", "", { "dependencies": { "reusify": "^1.0.4" } }, "sha512-YpkpUnK8od0o1hmeSc7UUs/eB/vIPWJYjKck2QKIzAf71Vm1AAQ3EbuZB3g2JIy+pg+ERD0vqI79KyZiB2e2Nw=="],
+
+    "fill-range": ["fill-range@7.0.1", "", { "dependencies": { "to-regex-range": "^5.0.1" } }, "sha512-qOo9F+dMUmC2Lcb4BbVvnKJxTPjCm+RRpe4gDuGrzkL7mEVl/djYSu2OdQ2Pa302N4oqkSg9ir6jaLWJ2USVpQ=="],
+
+    "filter-obj": ["filter-obj@1.1.0", "", {}, "sha1-mzERErxsYSehbgFsbF1/GeCAXFs="],
+
+    "form-data": ["form-data@4.0.0", "", { "dependencies": { "asynckit": "^0.4.0", "combined-stream": "^1.0.8", "mime-types": "^2.1.12" } }, "sha512-ETEklSGi5t0QMZuiXoA/Q6vcnxcLQP5vdugSpuAyi6SVGi2clPPp+xgEhuMaHC+zGgn31Kd235W35f7Hykkaww=="],
+
+    "fs-extra": ["fs-extra@9.1.0", "", { "dependencies": { "at-least-node": "^1.0.0", "graceful-fs": "^4.2.0", "jsonfile": "^6.0.1", "universalify": "^2.0.0" } }, "sha512-hcg3ZmepS30/7BSFqRvoo3DOMQu7IjqxO5nCDt+zM9XWjb33Wg7ziNT+Qvqbuc3+gWpzO02JubVyk2G4Zvo1OQ=="],
+
+    "fs.realpath": ["fs.realpath@1.0.0", "", {}, "sha1-FQStJSMVjKpA20onh8sBQRmU6k8="],
+
+    "function-bind": ["function-bind@1.1.1", "", {}, "sha512-yIovAzMX49sF8Yl58fSCWJ5svSLuaibPxXQJFLmBObTuCr0Mf1KiPopGM9NiFjiYBCbfaa2Fh6breQ6ANVTI0A=="],
+
+    "get-intrinsic": ["get-intrinsic@1.1.1", "", { "dependencies": { "function-bind": "^1.1.1", "has": "^1.0.3", "has-symbols": "^1.0.1" } }, "sha512-kWZrnVM42QCiEA2Ig1bG8zjoIMOgxWwYCEeNdwY6Tv/cOSeGpcoX4pXHfKUxNKVoArnrEr2e9srnAxxGIraS9Q=="],
+
+    "get-stream": ["get-stream@6.0.1", "", {}, "sha512-ts6Wi+2j3jQjqi70w5AlN8DFnkSwC+MqmxEzdEALB2qXZYV3X/b1CTfgPLGJNMeAWxdPfU8FO1ms3NUfaHCPYg=="],
+
+    "glob": ["glob@7.2.0", "", { "dependencies": { "fs.realpath": "^1.0.0", "inflight": "^1.0.4", "inherits": "2", "minimatch": "^3.0.4", "once": "^1.3.0", "path-is-absolute": "^1.0.0" } }, "sha512-lmLf6gtyrPq8tTjSmrO94wBeQbFR3HbLHbuyD69wuyQkImp2hWqMGB47OX65FBkPffO641IP9jWa1z4ivqG26Q=="],
+
+    "glob-parent": ["glob-parent@5.1.2", "", { "dependencies": { "is-glob": "^4.0.1" } }, "sha512-AOIgSQCepiJYwP3ARnGx+5VnTu2HBYdzbGP45eLw1vr3zB3vZLeyed1sC9hnbcOc9/SrMyM5RPQrkGz4aS9Zow=="],
+
+    "glob-promise": ["glob-promise@4.2.2", "", { "dependencies": { "@types/glob": "^7.1.3" }, "peerDependencies": { "glob": "^7.1.6" } }, "sha512-xcUzJ8NWN5bktoTIX7eOclO1Npxd/dyVqUJxlLIDasT4C7KZyqlPIwkdJ0Ypiy3p2ZKahTjK4M9uC3sNSfNMzw=="],
+
+    "globby": ["globby@11.1.0", "", { "dependencies": { "array-union": "^2.1.0", "dir-glob": "^3.0.1", "fast-glob": "^3.2.9", "ignore": "^5.2.0", "merge2": "^1.4.1", "slash": "^3.0.0" } }, "sha512-jhIXaOzy1sb8IyocaruWSn1TjmnBVs8Ayhcy83rmxNJ8q2uWKCAj3CnJY+KpGSXCueAPc0i05kVvVKtP1t9S3g=="],
+
+    "got": ["got@11.8.3", "", { "dependencies": { "@sindresorhus/is": "^4.0.0", "@szmarczak/http-timer": "^4.0.5", "@types/cacheable-request": "^6.0.1", "@types/responselike": "^1.0.0", "cacheable-lookup": "^5.0.3", "cacheable-request": "^7.0.2", "decompress-response": "^6.0.0", "http2-wrapper": "^1.0.0-beta.5.2", "lowercase-keys": "^2.0.0", "p-cancelable": "^2.0.0", "responselike": "^2.0.0" } }, "sha512-7gtQ5KiPh1RtGS9/Jbv1ofDpBFuq42gyfEib+ejaRBJuj/3tQFeR5+gw57e4ipaU8c/rCjvX6fkQz2lyDlGAOg=="],
+
+    "graceful-fs": ["graceful-fs@4.2.9", "", {}, "sha512-NtNxqUcXgpW2iMrfqSfR73Glt39K+BLwWsPs94yR63v45T0Wbej7eRmL5cWfwEgqXnmjQp3zaJTshdRW/qC2ZQ=="],
+
+    "has": ["has@1.0.3", "", { "dependencies": { "function-bind": "^1.1.1" } }, "sha512-f2dvO0VU6Oej7RkWJGrehjbzMAjFp5/VKPp5tTpWIV4JHHZK1/BxbFRtf/siA2SWTe09caDmVtYYzWEIbBS4zw=="],
+
+    "has-flag": ["has-flag@3.0.0", "", {}, "sha1-tdRU3CGZriJWmfNGfloH87lVuv0="],
+
+    "has-symbols": ["has-symbols@1.0.2", "", {}, "sha512-chXa79rL/UC2KlX17jo3vRGz0azaWEx5tGqZg5pO3NUyEJVB17dMruQlzCCOfUvElghKcm5194+BCRvi2Rv/Gw=="],
+
+    "http-cache-semantics": ["http-cache-semantics@4.1.0", "", {}, "sha512-carPklcUh7ROWRK7Cv27RPtdhYhUsela/ue5/jKzjegVvXDqM2ILE9Q2BGn9JZJh1g87cp56su/FgQSzcWS8cQ=="],
+
+    "http2-wrapper": ["http2-wrapper@1.0.3", "", { "dependencies": { "quick-lru": "^5.1.1", "resolve-alpn": "^1.0.0" } }, "sha512-V+23sDMr12Wnz7iTcDeJr3O6AIxlnvT/bmaAAAP/Xda35C90p9599p0F1eHR/N1KILWSoWVAiOMFjBBXaXSMxg=="],
+
+    "human-signals": ["human-signals@2.1.0", "", {}, "sha512-B4FFZ6q/T2jhhksgkbEW3HBvWIfDW85snkQgawt07S7J5QXTk6BkNV+0yAeZrM5QpMAdYlocGoljn0sJ/WQkFw=="],
+
+    "ignore": ["ignore@5.2.0", "", {}, "sha512-CmxgYGiEPCLhfLnpPp1MoRmifwEIOgjcHXxOBjv7mY96c+eWScsOP9c112ZyLdWHi0FxHjI+4uVhKYp/gcdRmQ=="],
+
+    "indent-string": ["indent-string@4.0.0", "", {}, "sha512-EdDDZu4A2OyIK7Lr/2zG+w5jmbuk1DVBnEwREQvBzspBJkCEbRa8GxU1lghYcaGJCnRWibjDXlq779X1/y5xwg=="],
+
+    "inflight": ["inflight@1.0.6", "", { "dependencies": { "once": "^1.3.0", "wrappy": "1" } }, "sha1-Sb1jMdfQLQwJvJEKEHW6gWW1bfk="],
+
+    "inherits": ["inherits@2.0.4", "", {}, "sha512-k/vGaX4/Yla3WzyMCvTQOXYeIHvqOKtnqBduzTHpzpQZzAskKMhZ2K+EnBiSM9zGSoIFeMpXKxa4dYeZIQqewQ=="],
+
+    "is-arrayish": ["is-arrayish@0.2.1", "", {}, "sha1-d8mYQFJ6qOyxqLppe4BkWnqSap0="],
+
+    "is-extglob": ["is-extglob@2.1.1", "", {}, "sha1-qIwCU1eR8C7TfHahueqXc8gz+MI="],
+
+    "is-glob": ["is-glob@4.0.3", "", { "dependencies": { "is-extglob": "^2.1.1" } }, "sha512-xelSayHH36ZgE7ZWhli7pW34hNbNl8Ojv5KVmkJD4hBdD3th8Tfk9vYasLM+mXWOZhFkgZfxhLSnrwRr4elSSg=="],
+
+    "is-number": ["is-number@7.0.0", "", {}, "sha512-41Cifkg6e8TylSpdtTpeLVMqvSBEVzTttHvERD741+pnZ8ANv0004MRL43QKPDlK9cGvNp6NZWZUBlbGXYxxng=="],
+
+    "is-ssh": ["is-ssh@1.3.3", "", { "dependencies": { "protocols": "^1.1.0" } }, "sha512-NKzJmQzJfEEma3w5cJNcUMxoXfDjz0Zj0eyCalHn2E6VOwlzjZo0yuO2fcBSf8zhFuVCL/82/r5gRcoi6aEPVQ=="],
+
+    "is-stream": ["is-stream@2.0.1", "", {}, "sha512-hFoiJiTl63nn+kstHGBtewWSKnQLpyb155KHheA1l39uvtO9nWIop1p3udqPcUd/xbF1VLMO4n7OI6p7RbngDg=="],
+
+    "isexe": ["isexe@2.0.0", "", {}, "sha1-6PvzdNxVb/iUehDcsFctYz8s+hA="],
+
+    "js-tokens": ["js-tokens@4.0.0", "", {}, "sha512-RdJUflcE3cUzKiMqQgsCu06FPu9UdIJO0beYbPhHN4k6apgJtifcoCtT9bcxOpYBtpD2kCM6Sbzg4CausW/PKQ=="],
+
+    "json-buffer": ["json-buffer@3.0.1", "", {}, "sha512-4bV5BfR2mqfQTJm+V5tPPdf+ZpuhiIvTuAB5g8kcrXOZpTT/QwwVRWBywX1ozr6lEuPdbHxwaJlm9G6mI2sfSQ=="],
+
+    "json-parse-even-better-errors": ["json-parse-even-better-errors@2.3.1", "", {}, "sha512-xyFwyhro/JEof6Ghe2iz2NcXoj2sloNsWr/XsERDK/oiPCfaNhl5ONfp+jQdAZRQQ0IJWNzH9zIZF7li91kh2w=="],
+
+    "jsonfile": ["jsonfile@6.1.0", "", { "dependencies": { "universalify": "^2.0.0" }, "optionalDependencies": { "graceful-fs": "^4.1.6" } }, "sha512-5dgndWOriYSm5cnYaJNhalLNDKOqFwyDB/rr1E9ZsGciGvKPs8R2xYGCacuf3z6K1YKDz182fd+fY3cn3pMqXQ=="],
+
+    "keyv": ["keyv@4.0.5", "", { "dependencies": { "json-buffer": "3.0.1" } }, "sha512-531pkGLqV3BMg0eDqqJFI0R1mkK1Nm5xIP2mM6keP5P8WfFtCkg2IOwplTUmlGoTgIg9yQYZ/kdihhz89XH3vA=="],
+
+    "lines-and-columns": ["lines-and-columns@1.2.4", "", {}, "sha512-7ylylesZQ/PV29jhEDl3Ufjo6ZX7gCqJr5F7PKrqc93v7fzSymt1BpwEU8nAUXs8qzzvqhbjhK5QZg6Mt/HkBg=="],
+
+    "lodash": ["lodash@4.17.21", "", {}, "sha512-v2kDEe57lecTulaDIuNTPy3Ry4gLGJ6Z1O3vE1krgXZNrsQ+LFTGHVxVjcXPs17LhbZVGedAJv8XZ1tvj5FvSg=="],
+
+    "lowercase-keys": ["lowercase-keys@2.0.0", "", {}, "sha512-tqNXrS78oMOE73NMxK4EMLQsQowWf8jKooH9g7xPavRT706R6bkQJ6DY2Te7QukaZsulxa30wQ7bk0pm4XiHmA=="],
+
+    "merge-stream": ["merge-stream@2.0.0", "", {}, "sha512-abv/qOcuPfk3URPfDzmZU1LKmuw8kT+0nIHvKrKgFrwifol/doWcdA4ZqsWQ8ENrFKkd67Mfpo/LovbIUsbt3w=="],
+
+    "merge2": ["merge2@1.4.1", "", {}, "sha512-8q7VEgMJW4J8tcfVPy8g09NcQwZdbwFEqhe/WZkoIzjn/3TGDwtOCYtXGxA3O8tPzpczCCDgv+P2P5y00ZJOOg=="],
+
+    "micromatch": ["micromatch@4.0.4", "", { "dependencies": { "braces": "^3.0.1", "picomatch": "^2.2.3" } }, "sha512-pRmzw/XUcwXGpD9aI9q/0XOwLNygjETJ8y0ao0wdqprrzDa4YnxLcz7fQRZr8voh8V10kGhABbNcHVk5wHgWwg=="],
+
+    "mime-db": ["mime-db@1.51.0", "", {}, "sha512-5y8A56jg7XVQx2mbv1lu49NR4dokRnhZYTtL+KGfaa27uq4pSTXkwQkFJl4pkRMyNFz/EtYDSkiiEHx3F7UN6g=="],
+
+    "mime-types": ["mime-types@2.1.34", "", { "dependencies": { "mime-db": "1.51.0" } }, "sha512-6cP692WwGIs9XXdOO4++N+7qjqv0rqxxVvJ3VHPh/Sc9mVZcQP+ZGhkKiTvWMQRr2tbHkJP/Yn7Y0npb3ZBs4A=="],
+
+    "mimic-fn": ["mimic-fn@2.1.0", "", {}, "sha512-OqbOk5oEQeAZ8WXWydlu9HJjz9WVdEIvamMCcXmuqUYjTknH/sqsWvhQ3vgwKFRR1HpjvNBKQ37nbJgYzGqGcg=="],
+
+    "mimic-response": ["mimic-response@3.1.0", "", {}, "sha512-z0yWI+4FDrrweS8Zmt4Ej5HdJmky15+L2e6Wgn3+iK5fWzb6T3fhNFq2+MeTRb064c6Wr4N/wv0DzQTjNzHNGQ=="],
+
+    "minimatch": ["minimatch@3.0.4", "", { "dependencies": { "brace-expansion": "^1.1.7" } }, "sha512-yJHVQEhyqPLUTgt9B83PXu6W3rx4MvvHvSUvToogpwoGDOUQ+yDrR0HRot+yOCdCO7u4hX3pWft6kWBBcqh0UA=="],
+
+    "ms": ["ms@2.1.2", "", {}, "sha512-sGkPx+VjMtmA6MX27oA4FBFELFCZZ4S4XqeGOXCv68tT+jb3vk/RyaKWP0PTKyWtmLSM0b+adUTEvbs1PEaH2w=="],
+
+    "normalize-url": ["normalize-url@6.1.0", "", {}, "sha512-DlL+XwOy3NxAQ8xuC0okPgK46iuVNAK01YN7RueYBqqFeGsBjV9XmCAzAdgt+667bCl5kPh9EqKKDwnaPG1I7A=="],
+
+    "npm-run-path": ["npm-run-path@4.0.1", "", { "dependencies": { "path-key": "^3.0.0" } }, "sha512-S48WzZW777zhNIrn7gxOlISNAqi9ZC/uQFnRdbeIHhZhCA6UqpkOT8T1G7BvfdgP4Er8gF4sUbaS0i7QvIfCWw=="],
+
+    "object-inspect": ["object-inspect@1.12.0", "", {}, "sha512-Ho2z80bVIvJloH+YzRmpZVQe87+qASmBUKZDWgx9cu+KDrX2ZDH/3tMy+gXbZETVGs2M8YdxObOh7XAtim9Y0g=="],
+
+    "once": ["once@1.4.0", "", { "dependencies": { "wrappy": "1" } }, "sha1-WDsap3WWHUsROsF9nFC6753Xa9E="],
+
+    "onetime": ["onetime@5.1.2", "", { "dependencies": { "mimic-fn": "^2.1.0" } }, "sha512-kbpaSSGJTWdAY5KPVeMOKXSrPtr8C8C7wodJbcsd51jRnmD+GZu8Y0VoU6Dm5Z4vWr0Ig/1NKuWRKf7j5aaYSg=="],
+
+    "p-cancelable": ["p-cancelable@2.1.1", "", {}, "sha512-BZOr3nRQHOntUjTrH8+Lh54smKHoHyur8We1V8DSMVrl5A2malOOwuJRnKRDjSnkoeBh4at6BwEnb5I7Jl31wg=="],
+
+    "p-reduce": ["p-reduce@2.1.0", "", {}, "sha512-2USApvnsutq8uoxZBGbbWM0JIYLiEMJ9RlaN7fAzVNb9OZN0SHjjTTfIcb667XynS5Y1VhwDJVDa72TnPzAYWw=="],
+
+    "parse-json": ["parse-json@5.2.0", "", { "dependencies": { "@babel/code-frame": "^7.0.0", "error-ex": "^1.3.1", "json-parse-even-better-errors": "^2.3.0", "lines-and-columns": "^1.1.6" } }, "sha512-ayCKvm/phCGxOkYRSCM82iDwct8/EonSEgCSxWxD7ve6jHggsFl4fZVQBPRNgQoKiuV/odhFrGzQXZwbifC8Rg=="],
+
+    "parse-path": ["parse-path@4.0.3", "", { "dependencies": { "is-ssh": "^1.3.0", "protocols": "^1.4.0", "qs": "^6.9.4", "query-string": "^6.13.8" } }, "sha512-9Cepbp2asKnWTJ9x2kpw6Fe8y9JDbqwahGCTvklzd/cEq5C5JC59x2Xb0Kx+x0QZ8bvNquGO8/BWP0cwBHzSAA=="],
+
+    "path-is-absolute": ["path-is-absolute@1.0.1", "", {}, "sha1-F0uSaHNVNP+8es5r9TpanhtcX18="],
+
+    "path-key": ["path-key@3.1.1", "", {}, "sha512-ojmeN0qd+y0jszEtoY48r0Peq5dwMEkIlCOu6Q5f41lfkswXuKtYrhgoTpLnyIcHm24Uhqx+5Tqm2InSwLhE6Q=="],
+
+    "path-type": ["path-type@4.0.0", "", {}, "sha512-gDKb8aZMDeD/tZWs9P6+q0J9Mwkdl6xMV8TjnGP3qJVJ06bdMgkbBlLU8IdfOsIsFz2BW1rNVT3XuNEl8zPAvw=="],
+
+    "picomatch": ["picomatch@2.3.1", "", {}, "sha512-JU3teHTNjmE2VCGFzuY8EXzCDVwEqB2a8fsIvwaStHhAWJEeVd1o1QD80CU6+ZdEXXSLbSsuLwJjkCBWqRQUVA=="],
+
+    "protocols": ["protocols@1.4.8", "", {}, "sha512-IgjKyaUSjsROSO8/D49Ab7hP8mJgTYcqApOqdPhLoPxAplXmkp+zRvsrSQjFn5by0rhm4VH0GAUELIPpx7B1yg=="],
+
+    "pump": ["pump@3.0.0", "", { "dependencies": { "end-of-stream": "^1.1.0", "once": "^1.3.1" } }, "sha512-LwZy+p3SFs1Pytd/jYct4wpv49HiYCqd9Rlc5ZVdk0V+8Yzv6jR5Blk3TRmPL1ft69TxP0IMZGJ+WPFU2BFhww=="],
+
+    "qs": ["qs@6.10.3", "", { "dependencies": { "side-channel": "^1.0.4" } }, "sha512-wr7M2E0OFRfIfJZjKGieI8lBKb7fRCH4Fv5KNPEs7gJ8jadvotdsS08PzOKR7opXhZ/Xkjtt3WF9g38drmyRqQ=="],
+
+    "query-string": ["query-string@6.14.1", "", { "dependencies": { "decode-uri-component": "^0.2.0", "filter-obj": "^1.1.0", "split-on-first": "^1.0.0", "strict-uri-encode": "^2.0.0" } }, "sha512-XDxAeVmpfu1/6IjyT/gXHOl+S0vQ9owggJ30hhWKdHAsNPOcasn5o9BW0eejZqL2e4vMjhAxoW3jVHcD6mbcYw=="],
+
+    "queue-microtask": ["queue-microtask@1.2.3", "", {}, "sha512-NuaNSa6flKT5JaSYQzJok04JzTL1CA6aGhv5rfLW3PgqA+M2ChpZQnAC8h8i4ZFkBS8X5RqkDBHA7r4hej3K9A=="],
+
+    "quick-lru": ["quick-lru@5.1.1", "", {}, "sha512-WuyALRjWPDGtt/wzJiadO5AXY+8hZ80hVpe6MyivgraREW751X3SbhRvG3eLKOYN+8VEvqLcf3wdnt44Z4S4SA=="],
+
+    "resolve-alpn": ["resolve-alpn@1.2.1", "", {}, "sha512-0a1F4l73/ZFZOakJnQ3FvkJ2+gSTQWz/r2KE5OdDY0TxPm5h4GkqkWWfM47T7HsbnOtcJVEF4epCVy6u7Q3K+g=="],
+
+    "responselike": ["responselike@2.0.0", "", { "dependencies": { "lowercase-keys": "^2.0.0" } }, "sha512-xH48u3FTB9VsZw7R+vvgaKeLKzT6jOogbQhEe/jewwnZgzPcnyWui2Av6JpoYZF/91uueC+lqhWqeURw5/qhCw=="],
+
+    "reusify": ["reusify@1.0.4", "", {}, "sha512-U9nH88a3fc/ekCF1l0/UP1IosiuIjyTh7hBvXVMHYgVcfGvt897Xguj2UOLDeI5BG2m7/uwyaLVT6fbtCwTyzw=="],
+
+    "run-parallel": ["run-parallel@1.2.0", "", { "dependencies": { "queue-microtask": "^1.2.2" } }, "sha512-5l4VyZR86LZ/lDxZTR6jqL8AFE2S0IFLMP26AbjsLVADxHdhB/c0GUsH+y39UfCi3dzz8OlQuPmnaJOMoDHQBA=="],
+
+    "sax": ["sax@1.2.4", "", {}, "sha512-NqVDv9TpANUjFm0N8uM5GxL36UgKi9/atZw+x7YFnQ8ckwFGKrl4xX4yWtrey3UJm5nP1kUbnYgLopqWNSRhWw=="],
+
+    "semantic-release-dotnet": ["semantic-release-dotnet@1.0.0", "", { "dependencies": { "glob": "^7.1.7", "glob-promise": "^4.2.0", "xml-js": "^1.6.11" } }, "sha512-U/cHwqqzFbJpPCQ/KMTSZtwzPNWCNCVStZRznMGK0xjhiLoDfRe5KFRs/9dzWBtPa358D7IiVzz97ZzepUAtfQ=="],
+
+    "shebang-command": ["shebang-command@2.0.0", "", { "dependencies": { "shebang-regex": "^3.0.0" } }, "sha512-kHxr2zZpYtdmrN1qDjrrX/Z1rR1kG8Dx+gkpK1G4eXmvXswmcE1hTWBWYUzlraYw1/yZp6YuDY77YtvbN0dmDA=="],
+
+    "shebang-regex": ["shebang-regex@3.0.0", "", {}, "sha512-7++dFhtcx3353uBaq8DDR4NuxBetBzC7ZQOhmTQInHEd6bSrXdiEyzCvG07Z44UYdLShWUyXt5M/yhz8ekcb1A=="],
+
+    "side-channel": ["side-channel@1.0.4", "", { "dependencies": { "call-bind": "^1.0.0", "get-intrinsic": "^1.0.2", "object-inspect": "^1.9.0" } }, "sha512-q5XPytqFEIKHkGdiMIrY10mvLRvnQh42/+GoBlFW3b2LXLE2xxJpZFdm94we0BaoV3RwJyGqg5wS7epxTv0Zvw=="],
+
+    "signal-exit": ["signal-exit@3.0.6", "", {}, "sha512-sDl4qMFpijcGw22U5w63KmD3cZJfBuFlVNbVMKje2keoKML7X2UzWbc4XrmEbDwg0NXJc3yv4/ox7b+JWb57kQ=="],
+
+    "slash": ["slash@3.0.0", "", {}, "sha512-g9Q1haeby36OSStwb4ntCGGGaKsaVSjQ68fBxoQcutl5fS1vuY18H3wSt3jFyFtrkx+Kz0V1G85A4MyAdDMi2Q=="],
+
+    "split-on-first": ["split-on-first@1.1.0", "", {}, "sha512-43ZssAJaMusuKWL8sKUBQXHWOpq8d6CfN/u1p4gUzfJkM05C8rxTmYrkIPTXapZpORA6LkkzcUulJ8FqA7Uudw=="],
+
+    "strict-uri-encode": ["strict-uri-encode@2.0.0", "", {}, "sha1-ucczDHBChi9rFC3CdLvMWGbONUY="],
+
+    "strip-final-newline": ["strip-final-newline@2.0.0", "", {}, "sha512-BrpvfNAE3dcvq7ll3xVumzjKjZQ5tI1sEUIKr3Uoks0XUl45St3FlatVqef9prk4jRDzhW6WZg+3bk93y6pLjA=="],
+
+    "supports-color": ["supports-color@5.5.0", "", { "dependencies": { "has-flag": "^3.0.0" } }, "sha512-QjVjwdXIt408MIiAqCX4oUKsgU2EqAGzs2Ppkm4aQYbjm+ZEWEcW4SfFNTr4uMNZma0ey4f5lgLrkB0aX0QMow=="],
+
+    "to-regex-range": ["to-regex-range@5.0.1", "", { "dependencies": { "is-number": "^7.0.0" } }, "sha512-65P7iz6X5yEr1cwcgvQxbbIw7Uk3gOy5dIdtZ4rDveLqhrdJP+Li/Hx6tyK0NEb+2GCyneCMJiGqrADCSNk8sQ=="],
+
+    "universalify": ["universalify@2.0.0", "", {}, "sha512-hAZsKq7Yy11Zu1DE0OzWjw7nnLZmJZYTDZZyEFHZdUhV8FkH5MCfoU1XMaxXovpyW5nq5scPqq0ZDP9Zyl04oQ=="],
+
+    "url-join": ["url-join@4.0.1", "", {}, "sha512-jk1+QP6ZJqyOiuEI9AEWQfju/nB2Pw466kbA0LEZljHwKeMgd9WrAEgEGxjPDD2+TNbbb37rTyhEfrCXfuKXnA=="],
+
+    "which": ["which@2.0.2", "", { "dependencies": { "isexe": "^2.0.0" }, "bin": { "node-which": "bin/node-which" } }, "sha512-BLI3Tl1TW3Pvl70l3yq3Y64i+awpwXqsGBYWkkqMtnbXgrMD+yj7rhW0kuEDxzJaYXGjEW5ogapKNMEKNMjibA=="],
+
+    "wrappy": ["wrappy@1.0.2", "", {}, "sha1-tSQ9jz7BqjXxNkYFvA0QNuMKtp8="],
+
+    "xml-js": ["xml-js@1.6.11", "", { "dependencies": { "sax": "^1.2.4" }, "bin": "bin/cli.js" }, "sha512-7rVi2KMfwfWFl+GpPg6m80IVMWXLRjO+PxTq7V2CDhoGak0wzYzFgUY2m4XJ47OGdXd8eLE8EmwfAmdjw7lC1g=="],
+
+    "@semantic-release/gitlab/fs-extra": ["fs-extra@10.0.0", "", { "dependencies": { "graceful-fs": "^4.2.0", "jsonfile": "^6.0.1", "universalify": "^2.0.0" } }, "sha512-C5owb14u9eJwizKGdchcDUQeFtlSHHthBk8pbX9Vc1PFZrLombudjDnNns88aYslCyF6IY5SUw3Roz6xShcEIQ=="],
+
+    "cacheable-request/get-stream": ["get-stream@5.2.0", "", { "dependencies": { "pump": "^3.0.0" } }, "sha512-nBF+F1rAZVCu/p7rjzgA+Yb4lfYXrpl7a6VmJrU8wF9I1CKvP/QwPNZHnOlwbTkY6dvtFIzFMSyQXbLoTQPRpA=="],
+
+    "chalk/escape-string-regexp": ["escape-string-regexp@1.0.5", "", {}, "sha1-G2HAViGQqN/2rjuyzwIAyhMLhtQ="],
+
+    "clone-response/mimic-response": ["mimic-response@1.0.1", "", {}, "sha512-j5EctnkH7amfV/q5Hgmoal1g2QHFJRraOtmx0JpIqkxhBhI/lJSl1nMpQ45hVarwNETOoWEimndZ4QK0RHxuxQ=="],
+  }
+}
--- a/default.nix
+++ b/default.nix
@@ -0,0 +1,19 @@
+{
+  sources ? import ./nix,
+  system ? builtins.currentSystem,
+  pkgs ? import sources.nixpkgs {
+    inherit system;
+    config = { };
+    overlays = [ ];
+  },
+}:
+let
+  sdk = pkgs.dotnetCorePackages.sdk_9_0;
+  sdslite = pkgs.callPackage ./nix/sdslite.nix { dotnet-sdk = sdk; };
+  projnetFsharp = pkgs.callPackage ./nix/projnet.fsharp.nix { dotnet-sdk = sdk; };
+in
+pkgs.callPackage ./src {
+  SDSLite = sdslite;
+  projnet = projnetFsharp;
+  dotnet-sdk = sdk;
+}
--- a/nix/default.nix
+++ b/nix/default.nix
@@ -0,0 +1,146 @@
+/*
+  This file is provided under the MIT licence:
+
+  Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+# Generated by npins. Do not modify; will be overwritten regularly
+let
+  data = builtins.fromJSON (builtins.readFile ./sources.json);
+  version = data.version;
+
+  # https://github.com/NixOS/nixpkgs/blob/0258808f5744ca980b9a1f24fe0b1e6f0fecee9c/lib/lists.nix#L295
+  range =
+    first: last: if first > last then [ ] else builtins.genList (n: first + n) (last - first + 1);
+
+  # https://github.com/NixOS/nixpkgs/blob/0258808f5744ca980b9a1f24fe0b1e6f0fecee9c/lib/strings.nix#L257
+  stringToCharacters = s: map (p: builtins.substring p 1 s) (range 0 (builtins.stringLength s - 1));
+
+  # https://github.com/NixOS/nixpkgs/blob/0258808f5744ca980b9a1f24fe0b1e6f0fecee9c/lib/strings.nix#L269
+  stringAsChars = f: s: concatStrings (map f (stringToCharacters s));
+  concatMapStrings = f: list: concatStrings (map f list);
+  concatStrings = builtins.concatStringsSep "";
+
+  # If the environment variable NPINS_OVERRIDE_${name} is set, then use
+  # the path directly as opposed to the fetched source.
+  # (Taken from Niv for compatibility)
+  mayOverride =
+    name: path:
+    let
+      envVarName = "NPINS_OVERRIDE_${saneName}";
+      saneName = stringAsChars (c: if (builtins.match "[a-zA-Z0-9]" c) == null then "_" else c) name;
+      ersatz = builtins.getEnv envVarName;
+    in
+    if ersatz == "" then
+      path
+    else
+      # this turns the string into an actual Nix path (for both absolute and
+      # relative paths)
+      builtins.trace "Overriding path of \"${name}\" with \"${ersatz}\" due to set \"${envVarName}\"" (
+        if builtins.substring 0 1 ersatz == "/" then
+          /. + ersatz
+        else
+          /. + builtins.getEnv "PWD" + "/${ersatz}"
+      );
+
+  mkSource =
+    name: spec:
+    assert spec ? type;
+    let
+      path =
+        if spec.type == "Git" then
+          mkGitSource spec
+        else if spec.type == "GitRelease" then
+          mkGitSource spec
+        else if spec.type == "PyPi" then
+          mkPyPiSource spec
+        else if spec.type == "Channel" then
+          mkChannelSource spec
+        else if spec.type == "Tarball" then
+          mkTarballSource spec
+        else
+          builtins.throw "Unknown source type ${spec.type}";
+    in
+    spec // { outPath = mayOverride name path; };
+
+  mkGitSource =
+    {
+      repository,
+      revision,
+      url ? null,
+      submodules,
+      hash,
+      branch ? null,
+      ...
+    }:
+    assert repository ? type;
+    # At the moment, either it is a plain git repository (which has an url), or it is a GitHub/GitLab repository
+    # In the latter case, there we will always be an url to the tarball
+    if url != null && !submodules then
+      builtins.fetchTarball {
+        inherit url;
+        sha256 = hash; # FIXME: check nix version & use SRI hashes
+      }
+    else
+      let
+        url =
+          if repository.type == "Git" then
+            repository.url
+          else if repository.type == "GitHub" then
+            "https://github.com/${repository.owner}/${repository.repo}.git"
+          else if repository.type == "GitLab" then
+            "${repository.server}/${repository.repo_path}.git"
+          else
+            throw "Unrecognized repository type ${repository.type}";
+        urlToName =
+          url: rev:
+          let
+            matched = builtins.match "^.*/([^/]*)(\\.git)?$" url;
+
+            short = builtins.substring 0 7 rev;
+
+            appendShort = if (builtins.match "[a-f0-9]*" rev) != null then "-${short}" else "";
+          in
+          "${if matched == null then "source" else builtins.head matched}${appendShort}";
+        name = urlToName url revision;
+      in
+      builtins.fetchGit {
+        rev = revision;
+        inherit name;
+        # hash = hash;
+        inherit url submodules;
+      };
+
+  mkPyPiSource =
+    { url, hash, ... }:
+    builtins.fetchurl {
+      inherit url;
+      sha256 = hash;
+    };
+
+  mkChannelSource =
+    { url, hash, ... }:
+    builtins.fetchTarball {
+      inherit url;
+      sha256 = hash;
+    };
+
+  mkTarballSource =
+    {
+      url,
+      locked_url ? url,
+      hash,
+      ...
+    }:
+    builtins.fetchTarball {
+      url = locked_url;
+      sha256 = hash;
+    };
+in
+if version == 5 then
+  builtins.mapAttrs mkSource data.pins
+else
+  throw "Unsupported format version ${toString version} in sources.json. Try running `npins upgrade`"
--- a/nix/projnet.fsharp.nix
+++ b/nix/projnet.fsharp.nix
@@ -0,0 +1,29 @@
+{
+  dotnet-sdk,
+  fetchFromGitLab,
+  buildDotnetModule
+}:
+let
+  src = fetchFromGitLab {
+    owner = "oceanbox";
+    repo = "ProjNet.FSharp";
+    # tag = "v5.2.0";
+    rev = "722ce0c23fda6844a81e995afbb2d81cbd5f38ec";
+    private = true;
+    forceFetchGit = true;
+    hash = "sha256-Rvnnf/D2x90pwgvTbXz307MJVBlVPK/cCf1hqj2VosE=";
+  };
+in
+buildDotnetModule {
+  name = "ProjNet.FSharp";
+
+  src = src;
+
+  dotnet-sdk = dotnet-sdk;
+
+  projectFile = "src/ProjNet.FSharp.fsproj";
+  nugetDeps = "${src}/src/deps.json";
+  packNupkg = true;
+
+  executables = [ ];
+}
--- a/nix/sdslite.nix
+++ b/nix/sdslite.nix
@@ -0,0 +1,29 @@
+{
+  dotnet-sdk,
+  fetchFromGitLab,
+  buildDotnetModule
+}:
+let
+  src = fetchFromGitLab {
+    owner = "oceanbox";
+    repo = "SDSlite";
+    # tag = "v2.8.0";
+    rev = "8c1a158206c37bc57a5bd726a792bd6a9cd2ec01";
+    private = true;
+    forceFetchGit = true;
+    hash = "sha256-i9pNrmH/VC0Q9FCldbWGdZHkqSL1cdYtAOs7vX+DlXM=";
+  };
+in
+buildDotnetModule {
+  name = "Oceanbox.SDSLite";
+
+  src = src;
+
+  dotnet-sdk = dotnet-sdk;
+
+  projectFile = "ScientificDataSet/ScientificDataSet.csproj";
+  nugetDeps = "${src}/ScientificDataSet/deps.json";
+  packNupkg = true;
+
+  executables = [ ];
+}
--- a/nix/sources.json
+++ b/nix/sources.json
@@ -0,0 +1,23 @@
+{
+  "pins": {
+    "nix-utils": {
+      "type": "Git",
+      "repository": {
+        "type": "Git",
+        "url": "https://git.sr.ht/~mrtz/nix-utils"
+      },
+      "branch": "trunk",
+      "submodules": false,
+      "revision": "098f594425d2b9dde0657becad0f6498d074f8b3",
+      "url": null,
+      "hash": "0hh52w1fkpr1xx6j8cjm6g88j2352yv2ysqm1q51j59y6f583vyb"
+    },
+    "nixpkgs": {
+      "type": "Channel",
+      "name": "nixpkgs-unstable",
+      "url": "https://releases.nixos.org/nixpkgs/nixpkgs-26.05pre905319.f720de590661/nixexprs.tar.xz",
+      "hash": "07n4hhch0j6n69b0zchdjg0l80z2xrdk7k57ykv90cvhklim5dz1"
+    }
+  },
+  "version": 5
+}
--- a/nix/sources.nix
+++ b/nix/sources.nix
@@ -0,0 +1,4 @@
+{
+    "ProjNet.FSharp" = "https://gitlab.com/api/v4/projects/35009572/packages/nuget/download";
+    "Oceanbox.SDSLite" = "https://gitlab.com/api/v4/projects/34025102/packages/nuget/download";
+}
--- a/package-lock.json
+++ b/package-lock.json
--- a/shell.nix
+++ b/shell.nix
@@ -0,0 +1,25 @@
+let
+  sources = import ./nix;
+  pkgs = import sources.nixpkgs {};
+  dotnet-sdk = pkgs.dotnetCorePackages.sdk_9_0;
+in
+with import <nixpkgs> { };
+mkShell rec {
+  packages = [
+    bun
+    dotnet-sdk
+    fsautocomplete
+    fantomas
+    npins
+    nixfmt
+    nuget-to-json
+  ];
+
+  buildInputs = [
+    netcdf
+    stdenv.cc.cc.lib
+  ];
+
+  DOTNET_ROOT = "${dotnet-sdk}/share/dotnet";
+  LD_LIBRARY_PATH = pkgs.lib.makeLibraryPath buildInputs;
+}
--- a/src/Adjoin.fs
+++ b/src/Adjoin.fs
@@ -1,7 +1,7 @@
 module Oceanbox.FvcomKit.Adjoin

 open FSharpPlus
-open FsKDTree
+//open FsKDTree
 open ProjNet.FSharp
 open Grid
 open Types
@@ -9,10 +9,10 @@ open Polygon
 open KdTree // C# version

 type CropBox = (float * float) * (float * float)
-type Mask = bool [,]
+type Mask = bool[,]
 type Pos = float * float
-type PosVec = Pos []
-type BiPos = float [,] * float [,]
+type PosVec = Pos[]
+type BiPos = float[,] * float[,]

 let private reProject (proj: IProj) ((lng, lat): BiPos) =
    lng
@@ -20,7 +20,7 @@ let private reProject (proj: IProj) ((lng, lat): BiPos) =
    |> Array2D.map proj.project
    |> unzip2D

-let private genCropIdx ((min, max): float * float) (m: float [,]) =
+let private genCropIdx ((min, max): float * float) (m: float[,]) =
    [|
        for i = 0 to Array2D.length1 m - 1 do
            for j = 0 to Array2D.length2 m - 1 do
@@ -44,22 +44,21 @@ let private mkCropBox (box: BBox) : CropBox =

 // make an array of indeces of wet roms grid points inside the fvcom domain
 let private genCullIdx (box: CropBox) (wet: Mask) (p: BiPos) =
-    genCropMask box p
-    |> Array.filter (fun (i, j) -> wet[i, j])
+    genCropMask box p |> Array.filter (fun (i, j) -> wet[i, j])

 // crop roms grid coordinates based in index mask of active, overlapping points
-let private cullCoords (p: BiPos) (idx: (int * int) []) =
+let private cullCoords (p: BiPos) (idx: (int * int)[]) =
    let x, y = p
    let x' = idx |> Array.map (fun (i, j) -> x[i, j])
    let y' = idx |> Array.map (fun (i, j) -> y[i, j])
    Array.zip x' y'

-let private createTree (points: Leaf<float, int>[]) =
-    let tree = KdTree<float, int>(2, KdTree.Math.DoubleMath())
+let private createTree (points: Leaf<int>[]) =
+    let tree = KdTree<float, int> (2, KdTree.Math.DoubleMath ())
    points
-    |> Array.iter (fun a -> tree.Add([| a.Pos.X; a.Pos.Y |], a.Data) |> ignore)
+    |> Array.iter (fun a -> tree.Add ([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
    if points.Length > 0 then
-        tree.Balance()
+        tree.Balance ()
    else
        // Log.Warning $"Empty kd-tree"
        ()
@@ -68,14 +67,14 @@ let private createTree (points: Leaf<float, int>[]) =
 // make a kd-tree for fast nearest neighbour lookup
 let private buildTree (points: (float * float) array) =
    points
-    |> Array.mapi (fun n (x, y) -> { Pos = { X = x; Y = y }; Data = n })
+    |> Array.mapi (fun n (x, y) -> { Pos = x, y; Data = n })
    // |> create2DTree treeLeafSize
    |> createTree

 // adjoinIdx: index to culled roms grid for each fvcom node
 // cullIdx: index to roms points _actually_ in use
 // oobIdx: out-of-bounds points, too far from a roms cell
-type FvcomAdjoint = { adjoinIdx: int []; cullIdx: (int * int) []; oobIdx: int [] }
+type FvcomAdjoint = { adjoinIdx: int[]; cullIdx: (int * int)[]; oobIdx: int[] }

 let private dist (a: float * float) (b: float * float) =
    let x0, y0 = a
@@ -107,14 +106,13 @@ let private distLngLat (proj: IProj) (a: float * float) (b: float * float) =
 //         (0, [])
 //     |> snd
 //     |> Array.ofList
-let private genOobIdx (proj: IProj) (tree: KdTree<_, _>) (cullIdx: (int * int) []) (fPos: PosVec) (rPos: BiPos) =
+let private genOobIdx (proj: IProj) (tree: KdTree<_, _>) (cullIdx: (int * int)[]) (fPos: PosVec) (rPos: BiPos) =
    let rLon, rLat = rPos
-    let dMax =
-        distLngLat proj (rLon[0, 0], rLat[0, 0]) (rLon[0, 1], rLat[0, 1])
+    let dMax = distLngLat proj (rLon[0, 0], rLat[0, 0]) (rLon[0, 1], rLat[0, 1])
    fPos
    |> Array.fold
        (fun (n, a) (x, y as p0) ->
-            tree.GetNearestNeighbours([|x; y|], 1)
+            tree.GetNearestNeighbours ([| x; y |], 1)
            |> Array.head
            |> fun p ->
                let i0, i1 = cullIdx[p.Value]
@@ -141,6 +139,7 @@ let private genOobIdx (proj: IProj) (tree: KdTree<_, _>) (cullIdx: (int * int) [
 //         cullIdx = cullIdx
 //         oobIdx = genOobIdx proj tree cullIdx fPos rPos
 //     }
+
 let private mkFvcomAdjoint (proj: IProj) (fPos: PosVec) (bbox: BBox) ((rPos, wet): BiPos * Mask) =
    let box = mkCropBox bbox
    let pos' = reProject proj rPos
@@ -148,15 +147,10 @@ let private mkFvcomAdjoint (proj: IProj) (fPos: PosVec) (bbox: BBox) ((rPos, wet
    let tree = cullCoords pos' cullIdx |> buildTree
    let nearest =
        fPos
-        |> Array.map (fun (x, y) ->
-            tree.GetNearestNeighbours([|x; y|], 1)
-            |> Array.head
-            |> fun x -> x.Value)
-    {
-        adjoinIdx = nearest
-        cullIdx = cullIdx
-        oobIdx = genOobIdx proj tree cullIdx fPos rPos
-    }
+        |> Array.map (fun (x, y) -> tree.GetNearestNeighbours ([| x; y |], 1) |> Array.head |> (fun x -> x.Value))
+    let oobIdx = genOobIdx proj tree cullIdx fPos rPos
+
+    { adjoinIdx = nearest; cullIdx = cullIdx; oobIdx = oobIdx }

 let inline float2 x = bimap float float x

@@ -169,67 +163,58 @@ let getNearestCell (proj: IProj) (fvcom: Grid) (roms: BiPos * Mask) =
    let cells = Util.calcCentroids fvcom |> Array.map float2
    mkFvcomAdjoint proj cells fvcom.BBox roms

-let private createIdxTree (points: Leaf<float, int * int>[]) =
-    let tree = KdTree<float, int * int>(2, KdTree.Math.DoubleMath())
+let private createIdxTree (points: Leaf<int * int>[]) =
+    let tree = KdTree<float, int * int> (2, KdTree.Math.DoubleMath ())
    points
-    |> Array.iter (fun a -> tree.Add([| a.Pos.X; a.Pos.Y |], a.Data) |> ignore)
+    |> Array.iter (fun a -> tree.Add ([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
    if points.Length > 0 then
-        tree.Balance()
+        tree.Balance ()
    else
        //Log.Warning $"Empty kd-tree"
        ()
    tree

 // rPos must be in the right projection!
-let makeNestTree ((lng, lat): float [,] * float [,]) =
+let makeNestTree ((lng, lat): float[,] * float[,]) =
    [| 0 .. (Array2D.length1 lng) - 1 |]
    |> Array.map (fun i ->
        [| 0 .. (Array2D.length2 lng) - 1 |]
        |> Array.map (fun j ->
            let x = lng[i, j]
            let y = lat[i, j]
-            { Pos = { X = x; Y = y }; Data = (i, j) }))
+            { Pos = x, y; Data = (i, j) }))
    |> Array.concat
    // |> create2DTree treeLeafSize
    |> createIdxTree

-let getNearestUpperLeft'
-    (tree: KdTree<float, int * int>)
-    (((ew, ns), mask): BiPos * Mask)
-    ((x, y): float * float)
-    =
+let getNearestUpperLeft' (tree: KdTree<float, int * int>) (((ew, ns), mask): BiPos * Mask) ((x, y): float * float) =
    let i, j =
        // nearestNeighbor tree { X = x; Y = y }
-        tree.GetNearestNeighbours([|x; y|], 1)
-        |> fun x -> x[0].Value
-    let p1 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i + 1, j], ns[i + 1, j]
-            ew[i + 1, j - 1], ns[i + 1, j - 1]
-            ew[i, j - 1], ns[i, j - 1]
-        |] // lower right grid cell
-    let p2 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i + 1, j], ns[i + 1, j]
-            ew[i + 1, j + 1], ns[i + 1, j + 1]
-            ew[i, j + 1], ns[i, j + 1]
-        |] // upper right grid cell
-    let p3 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i - 1, j], ns[i - 1, j]
-            ew[i - 1, j - 1], ns[i - 1, j - 1]
-            ew[i, j - 1], ns[i, j - 1]
-        |] // lower left grid cell
-    let p4 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i - 1, j], ns[i - 1, j]
-            ew[i - 1, j + 1], ns[i - 1, j + 1]
-            ew[i, j + 1], ns[i, j + 1]
-        |] // upper left grid cell
+        tree.GetNearestNeighbours ([| x; y |], 1) |> fun x -> x[0].Value
+    let p1 = [|
+        ew[i, j], ns[i, j]
+        ew[i + 1, j], ns[i + 1, j]
+        ew[i + 1, j - 1], ns[i + 1, j - 1]
+        ew[i, j - 1], ns[i, j - 1]
+    |] // lower right grid cell
+    let p2 = [|
+        ew[i, j], ns[i, j]
+        ew[i + 1, j], ns[i + 1, j]
+        ew[i + 1, j + 1], ns[i + 1, j + 1]
+        ew[i, j + 1], ns[i, j + 1]
+    |] // upper right grid cell
+    let p3 = [|
+        ew[i, j], ns[i, j]
+        ew[i - 1, j], ns[i - 1, j]
+        ew[i - 1, j - 1], ns[i - 1, j - 1]
+        ew[i, j - 1], ns[i, j - 1]
+    |] // lower left grid cell
+    let p4 = [|
+        ew[i, j], ns[i, j]
+        ew[i - 1, j], ns[i - 1, j]
+        ew[i - 1, j + 1], ns[i - 1, j + 1]
+        ew[i, j + 1], ns[i, j + 1]
+    |] // upper left grid cell
    let q =
        if inpolygon p1 (x, y) then
            [| i + 1, j - 1; i + 1, j; i, j; i, j - 1 |]
@@ -237,85 +222,88 @@ let getNearestUpperLeft'
            [| i + 1, j; i + 1, j + 1; i, j + 1; i, j |]
        elif inpolygon p3 (x, y) then
            [| i, j - 1; i, j; i - 1, j; i - 1, j - 1 |]
-        elif inpolygon p4 (x,y) then
+        elif inpolygon p4 (x, y) then
            [| i, j; i, j + 1; i - 1, j + 1; i - 1, j |]
        else
-            [| 0,0; 0,0; 0,0; 0,0 |]
+            [| 0, 0; 0, 0; 0, 0; 0, 0 |]
    if (q |> Array.map (fun (q1, q2) -> q1 + q2) |> Array.sum) = 0 then
        failwith "Surrounding grid cell not found"
    let m = q |> Array.map (fun (n, m) -> mask[n, m])
    Array.zip q m

-let getNearestUpperLeft
-    (tree: Tree<Leaf<float, int * int> array, Node<float>>)
-    (((ew, ns), mask): BiPos * Mask)
-    ((x, y): float * float)
-    =
-    let i, j =
-        nearestNeighbor tree { X = x; Y = y }
-        |> fun x -> x.Value.Data
-    let p1 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i + 1, j], ns[i + 1, j]
-            ew[i + 1, j - 1], ns[i + 1, j - 1]
-            ew[i, j - 1], ns[i, j - 1]
-        |] // lower right grid cell
-    let p2 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i + 1, j], ns[i + 1, j]
-            ew[i + 1, j + 1], ns[i + 1, j + 1]
-            ew[i, j + 1], ns[i, j + 1]
-        |] // upper right grid cell
-    let p3 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i - 1, j], ns[i - 1, j]
-            ew[i - 1, j - 1], ns[i - 1, j - 1]
-            ew[i, j - 1], ns[i, j - 1]
-        |] // lower left grid cell
-    let p4 =
-        [|
-            ew[i, j], ns[i, j]
-            ew[i - 1, j], ns[i - 1, j]
-            ew[i - 1, j + 1], ns[i - 1, j + 1]
-            ew[i, j + 1], ns[i, j + 1]
-        |] // upper left grid cell
-    let q =
-        if inpolygon p1 (x, y) then
-            [| i + 1, j - 1; i + 1, j; i, j; i, j - 1 |]
-        elif inpolygon p2 (x, y) then
-            [| i + 1, j; i + 1, j + 1; i, j + 1; i, j |]
-        elif inpolygon p3 (x, y) then
-            [| i, j - 1; i, j; i - 1, j; i - 1, j - 1 |]
-        elif inpolygon p4 (x,y) then
-            [| i, j; i, j + 1; i - 1, j + 1; i - 1, j |]
-        else
-            [| 0,0; 0,0; 0,0; 0,0 |]
-    if (q |> Array.map (fun (q1, q2) -> q1 + q2) |> Array.sum) = 0 then
-        failwith "Surrounding grid cell not found"
-    let m = q |> Array.map (fun (n, m) -> mask[n, m])
-    Array.zip q m
+// TODO: this is the F# Kd-tree version, which is currently defunct
+// let getNearestUpperLeft
+//     (tree: Tree<Leaf<float, int * int> [], Node<float>>)
+//     (((ew, ns), mask): BiPos * Mask)
+//     ((x, y): float * float)
+//     =
+//     let i, j =
+//         nearestNeighbor tree { X = x; Y = y }
+//         |> fun x -> x.Value.Data
+//     let p1 =
+//         [|
+//             ew[i, j], ns[i, j]
+//             ew[i + 1, j], ns[i + 1, j]
+//             ew[i + 1, j - 1], ns[i + 1, j - 1]
+//             ew[i, j - 1], ns[i, j - 1]
+//         |] // lower right grid cell
+//     let p2 =
+//         [|
+//             ew[i, j], ns[i, j]
+//             ew[i + 1, j], ns[i + 1, j]
+//             ew[i + 1, j + 1], ns[i + 1, j + 1]
+//             ew[i, j + 1], ns[i, j + 1]
+//         |] // upper right grid cell
+//     let p3 =
+//         [|
+//             ew[i, j], ns[i, j]
+//             ew[i - 1, j], ns[i - 1, j]
+//             ew[i - 1, j - 1], ns[i - 1, j - 1]
+//             ew[i, j - 1], ns[i, j - 1]
+//         |] // lower left grid cell
+//     let p4 =
+//         [|
+//             ew[i, j], ns[i, j]
+//             ew[i - 1, j], ns[i - 1, j]
+//             ew[i - 1, j + 1], ns[i - 1, j + 1]
+//             ew[i, j + 1], ns[i, j + 1]
+//         |] // upper left grid cell
+//     let q =
+//         if inpolygon p1 (x, y) then
+//             [| i + 1, j - 1; i + 1, j; i, j; i, j - 1 |]
+//         elif inpolygon p2 (x, y) then
+//             [| i + 1, j; i + 1, j + 1; i, j + 1; i, j |]
+//         elif inpolygon p3 (x, y) then
+//             [| i, j - 1; i, j; i - 1, j; i - 1, j - 1 |]
+//         elif inpolygon p4 (x, y) then
+//             [| i, j; i, j + 1; i - 1, j + 1; i - 1, j |]
+//         else
+//             [| 0, 0; 0, 0; 0, 0; 0, 0 |]
+//     if (q
+//         |> Array.map (fun (q1, q2) -> q1 + q2)
+//         |> Array.sum) = 0 then
+//         failwith "Surrounding grid cell not found"
+//     let m = q |> Array.map (fun (n, m) -> mask[n, m])
+//     Array.zip q m

-let getNearestCellCorner (coords, _ as grid: BiPos * Mask) (pos: (float * float) []) =
+let getNearestCellCorner (coords, _ as grid: BiPos * Mask) (pos: (float * float)[]) =
    let tree = makeNestTree coords
    // pos |> Array.map (getNearestUpperLeft tree grid)
    pos |> Array.map (getNearestUpperLeft' tree grid)

-let getCellBox ((ew, ns): BiPos) (boxid: ((int * int) * bool) []) =
+let getCellBox ((ew, ns): BiPos) (boxid: ((int * int) * bool)[]) =
    boxid
    |> Array.map (fun ((n, m), mask) -> (ew[n, m], ns[n, m]), mask)
    |> Array.unzip

-let getCellProps (prop: float [,]) ((n, m): int * int) =
+let getCellProps (prop: float[,]) ((n, m): int * int) =
    prop[n, m], prop[n + 1, m], prop[n + 1, m + 1], prop[n + 1, m]

 // pick out elements actually in use
-let inline private cullBiMatrix (culler: (int * int) []) (m: 'a [,]) =
+let inline private cullBiMatrix (culler: (int * int)[]) (m: 'a[,]) =
    culler |> Array.map (fun (i, j) -> m[i, j])

 // adjoin fvcom and roms data based on nearest neighbours
-let adjoinMatrix (adj: FvcomAdjoint) (m: 'a [,]) =
+let adjoinMatrix (adj: FvcomAdjoint) (m: 'a[,]) =
    let x = cullBiMatrix adj.cullIdx m
-    adj.adjoinIdx |> Array.map (fun n -> x[n])
+    adj.adjoinIdx |> Array.map (fun n -> x[n])
--- a/src/Arome.fs
+++ b/src/Arome.fs
@@ -0,0 +1,260 @@
+module Oceanbox.FvcomKit.Arome
+
+open System
+
+open Microsoft.Research.Science.Data
+
+open ProjNet.FSharp
+open Serilog
+
+open Types
+
+let trans = makeTransform CoordSys.WGS84 (CoordSys.LCCMet ())
+
+[<Struct>]
+type Pointf = {
+    x: float
+    y: float
+} with
+    static member Zero = { x = 0.0; y = 0.0 }
+    static member OfTuple (x, y) = { x = x; y = y }
+    static member OfStructTuple struct (x, y) = { x = x; y = y }
+
+/// Single precision
+[<Struct>]
+type Points = {
+    x: single
+    y: single
+} with
+    static member Zero = { x = 0.0f; y = 0.0f }
+    static member OfTuple (x, y) = { x = x; y = y }
+    static member OfStructTuple struct (x, y) = { x = x; y = y }
+
+module Pointf =
+    let ofPoints (p: Points) : Pointf = { x = float p.x; y = float p.y }
+
+    let getBBox (points: Pointf array) : BBox =
+        let minX = points |> Array.minBy _.x |> _.x
+        let maxX = points |> Array.maxBy _.x |> _.x
+        let minY = points |> Array.minBy _.y |> _.y
+        let maxY = points |> Array.maxBy _.y |> _.y
+        let center = (minX + (maxX - minX)) / 2., (minY + (maxY - minY)) / 2.
+
+        {
+            minX = minX
+            maxX = maxX
+            minY = minY
+            maxY = maxY
+            center = center
+        }
+
+module Points =
+    let getBBox (points: Points array) : BBox =
+        let minX = points |> Array.minBy _.x |> _.x
+        let maxX = points |> Array.maxBy _.x |> _.x
+        let minY = points |> Array.minBy _.y |> _.y
+        let maxY = points |> Array.maxBy _.y |> _.y
+        let center = float (minX + (maxX - minX)) / 2., float (minY + (maxY - minY)) / 2.
+
+        {
+            minX = float minX
+            maxX = float maxX
+            minY = float minY
+            maxY = float maxY
+            center = center
+        }
+
+[<Struct>]
+type SquareGrid = {
+    dimensions: int * int
+    BBox: BBox
+    squareSize: float
+    projection: ProjNet.CoordinateSystems.CoordinateSystem
+    points: Points array
+} with
+    member this.getBoundingBox() = this.BBox
+    static member empty = {
+        dimensions = 0, 0
+        BBox = BBox.empty
+        squareSize = 0.0
+        projection = CoordSys.LCCMet ()
+        points = Array.empty
+    }
+
+let private square x = x * x
+
+let haversineDistance (earthRadius: float) (x0: float) (y0: float) (x1: float) (y1: float) : float =
+    let mutable lat1 = y0
+    let mutable lat2 = y1
+    let lon1 = x0
+    let lon2 = x1
+
+    let dLat = Double.DegreesToRadians (lat2 - lat1)
+    let dLon = Double.DegreesToRadians (lon2 - lon1)
+    lat1 <- Double.DegreesToRadians lat1
+    lat2 <- Double.DegreesToRadians lat2
+
+    let a = square (Math.Sin (dLat / 2.0)) + Math.Cos lat1 * Math.Cos lat2 * square (Math.Sin (dLon / 2.0))
+    let c = 2.0 * Math.Asin (Math.Sqrt a)
+
+    let result = earthRadius * c
+
+    result
+
+let getGrid (ds: DataSet) : Result<SquareGrid, string> =
+    try
+        let dimensions = ds.Dimensions["x"].Length, ds.Dimensions["y"].Length
+        let longs : float array2d = ds["longitude"].GetData () :?> float[,]
+        let lats : float array2d = ds["latitude"].GetData () :?> float[,]
+        // NOTE: The netcdf file dimensions are defined as (y, x)
+        let width = Array2D.length2 longs
+        let height = Array2D.length1 lats
+        let points : Points array =
+            let result = Array.create (width * height) Points.Zero
+            for i in 0 .. height - 1 do
+                for j in 0 .. width - 1 do
+                    let lat =  lats[i, j]
+                    let lon = longs[i, j]
+                    let p = lon, lat
+                    // NOTE(simkir): Convert to lambert projection
+                    let x, y = trans.project p
+                    result[i * width + j] <- Points.OfTuple (single x, single y)
+            result
+
+        let bbox = Points.getBBox points
+
+        if points.Length < 2 then
+            Error "The dataset must contain at least 1 square"
+        else
+            let p0 = points[0 * width + 0]
+            let p1 = points[0 * width + 1]
+            let p2 = points[1 * width + 0]
+            // let p3 = points[1 * width + 1]
+            let x1, x0 = if p1.x > p0.x then p1.x, p0.x else p0.x, p1.x
+            let y1, y0 = if p2.y > p1.y then p2.y, p1.y else p1.y, p2.y
+            let lengthX = x1 - x0
+            let lengthY = y1 - y0
+            let isSquare = lengthX = lengthY
+
+            if not isSquare then
+                Log.Warning (
+                    "FvcomKit.Arome.getGrid grid is not square: {X1} - {X0} = {LengthX} = {LengthY} = {Y1} - {Y0}",
+                    x1,
+                    x0,
+                    lengthX,
+                    lengthY,
+                    y1,
+                    y0
+                )
+
+            Ok {
+                SquareGrid.empty with
+                    dimensions = dimensions
+                    BBox = bbox
+                    squareSize = float lengthX
+                    points = points
+            }
+    with exn ->
+        Log.Error (exn, "Sorcerer.Arome.getAromeSquareGrid exception")
+        Error $"Error reading arome grid: {exn.Message}"
+
+let private getBBox (xs: float array) (ys: float array) : BBox =
+    try
+        let minX = xs |> Array.min
+        let maxX = xs |> Array.max
+        let minY = ys |> Array.min
+        let maxY = ys |> Array.max
+        let center = (minX + (maxX - minX)) / 2., (minY + (maxY - minY)) / 2.
+
+        {
+            minX = minX
+            maxX = maxX
+            minY = minY
+            maxY = maxY
+            center = center
+        }
+    with e ->
+        Log.Error $"{e}"
+        BBox.empty
+
+/// Depends on the netcdf having the 'x' and 'y' variables
+let getSquareGrid (ds: DataSet) : Result<SquareGrid, string> =
+    try
+        let dimensions = ds.Dimensions["x"].Length, ds.Dimensions["y"].Length
+        let xs = (ds["x"].GetData () :?> single[]) |> Array.map float
+        let ys = (ds["y"].GetData () :?> single[]) |> Array.map float
+        let bbox = getBBox xs ys
+
+        if xs.Length < 2 && ys.Length < 2 then
+            Error "The dataset must contain at least 1 square"
+        else
+            let x0, y0 = xs[0], ys[0]
+            let x1, y1 = xs[1], ys[1]
+            let lengthX = if x1 > x0 then x1 - x0 else x0 - x1
+            let lengthY = if y1 > y0 then y1 - y0 else y0 - y1
+            let isSquare = lengthX = lengthY
+
+            if isSquare then
+                { SquareGrid.empty with dimensions = dimensions; BBox = bbox; squareSize = lengthX }
+                |> Ok
+            else
+                Log.Error (
+                    "FvcomKit.Arome.getGrid grid is not square: {X1} - {X0} = {LengthX} = {LengthY} = {Y1} - {Y0}",
+                    x1,
+                    x0,
+                    lengthX,
+                    lengthY,
+                    y1,
+                    y0
+                )
+                Error "The given data set does not contain a grid made up of squares"
+    with exn ->
+        Log.Error (exn, "FvcomKit.Arome.getGrid exception")
+        Error $"Error reading arome grid: {exn.Message}"
+
+let readUV (ds: DataSet) (t: int) x y : single * single =
+    let xWind = ds["x_wind_10m"]
+    let yWind = ds["y_wind_10m"]
+    let origin = [| t; 0; y; x |]
+    let shape = [| 1; 1; 1; 1; |]
+    Log.Verbose("""Fetching NetCDF["x_wind_10m"]({Origin}, {Shape})""", origin, shape)
+    let us = xWind.GetData (origin, shape) :?> single[,,,]
+    let vs = yWind.GetData (origin, shape) :?> single[,,,]
+
+    us[0, 0, 0, 0], vs[0, 0, 0, 0]
+
+/// Finds the index of a tile within a square grid, given its bounding box and square length
+let tryFindIndex (grid: SquareGrid) (x0, y0) =
+    let wide, tall = grid.dimensions
+    let bbox = grid.BBox
+    if bbox.minX <= x0 && x0 < bbox.maxX && bbox.minY <= y0 && y0 < bbox.maxY then
+        let dx = x0 - bbox.minX
+        let dy = y0 - bbox.minY
+        let xIdx = int (dx / grid.squareSize)
+        let yIdx = int (dy / grid.squareSize)
+
+        if xIdx < wide && yIdx < tall then
+            Some (xIdx, yIdx)
+        else
+            Log.Warning (
+                "Got wrong indices within the bounding box of the archive: min {@Min}, max {@Max}, point {@Point}, delta {@Delta}m, indices {@Indices}",
+                (bbox.minX, bbox.minY),
+                (bbox.maxX, bbox.maxY),
+                (x0, y0),
+                (dx, dy),
+                (xIdx, yIdx)
+            )
+            None
+    else
+        None
+
+/// Tries to get the closest x and y in the arome dataset based on position p
+let tryFind (grid: SquareGrid) (p: float * float) : (int * int) option =
+    tryFindIndex grid p
+
+let tryFindWithProj (proj: Projection) (grid: SquareGrid) (p0: float, p1: float) : (int * int) option =
+    let coordSys : ProjNet.CoordinateSystems.CoordinateSystem = Projection.ToCoordinateSystem proj
+    let trans = makeTransform coordSys grid.projection
+    let p = trans.project ((p0, p1))
+
+    tryFindIndex grid p
--- a/src/Evaluate.fs
+++ b/src/Evaluate.fs
@@ -1,236 +0,0 @@
-module Oceanbox.FvcomKit.Evaluate
-
-open Grid
-open Gradient
-
-let private getCell (grid: ExtendedGrid) (idx: ElemIdx) =
-    (grid :> IGrid).getCell idx
-
-let private getVertex (grid: ExtendedGrid) (idx: NodeIdx) =
-    (grid :> IGrid).getVertex idx
-
-let private getCentroid (grid: ExtendedGrid) (idx: ElemIdx) =
-    (grid.getCentroids ())[idx]
-
-let private containsIndex (idx: ElemIdx) (cell: Cell) =
-    let x, y, z = cell
-    (idx = x) || (idx = y) || (idx = z)
-
-let private countIndex (idx: ElemIdx) (cell: Cell) =
-    match (containsIndex idx cell) with
-    | true -> 1
-    | false -> 0
-
-let private countCommonIndices (c1: Cell) (c2: Cell) =
-    let i0, i1, i2 = c1
-    let s0 = countIndex i0 c2
-    let s1 = countIndex i1 c2
-    let s2 = countIndex i2 c2
-    s0 + s1 + s2
-
-// Elements are adjacent if they have a shared edge, i.e. have two common node indices
-let private areAdjacent (grid: ExtendedGrid) (i1: ElemIdx) (i0: ElemIdx) =
-    let nCommon =
-        countCommonIndices
-        <| (getCell grid i0)
-        <| (getCell grid i1)
-    nCommon = 2
-
-// Picks out the three neighbors of an element that are sharing a common edge
-// (two neighbors for boundary elements)
-let private getAdjacentNeighbors (grid: ExtendedGrid) (idx: ElemIdx) =
-    let neighbors =
-        grid.getElemsSurroundingElem idx
-        |> Array.filter (areAdjacent grid idx)
-
-    // a valid centroid grid should always give two or three neighbors
-    match neighbors.Length with
-    | n when n < 2 -> failwith "Too few neighbors"
-    | n when n > 3 -> failwith "Too many neighbors"
-    | _ -> neighbors
-
-// Evaluate 2D field based on linear extrapolation from centroid field value,
-// using gradient estimated from values in the three adjacent neigboring elements
-let evaluateLinearCentroid (grid: ExtendedGrid) (readUV: int -> single * single) (p: Pos) e0 (u, v) =
-    let centroids = grid.getCentroids ()
-
-    let c_0 = centroids[e0]
-    let u_0 = u
-    let v_0 = v
-
-    let neighbors = getAdjacentNeighbors grid e0
-    let c_n = neighbors |> Array.map (Array.get centroids)
-    let u_n = neighbors |> Array.map (fun eIdx -> readUV eIdx |> fst)
-    let v_n = neighbors |> Array.map (fun eIdx -> readUV eIdx |> snd)
-
-    let thrs = 1e-6f
-    let du_x, du_y = calcGradientFromNeighbors c_0 u_0 c_n u_n thrs
-    let dv_x, dv_y = calcGradientFromNeighbors c_0 v_0 c_n v_n thrs
-
-    let x, y = p
-    let x_0, y_0 = c_0
-
-    let u = u_0 + du_x*(x - x_0) + du_y*(y - y_0)
-    let v = v_0 + dv_x*(x - x_0) + dv_y*(y - y_0)
-    Some (u, v)
-
-// Evaluate 2D field based on constant extrapolation from centroid field value
-let evaluateConstantCentroid (grid: ExtendedGrid) (field: Field) (p: Pos)  =
-    match grid.tryGetElement p with
-        | None -> None
-        | Some e0 -> Some (field[e0])
-
-
-
-
-
-// Given centroid 2D field values, compute 2D field value in single vertex
-// by a weighted average of values in all surrounding elements
-let private calcVertexValue (grid: ExtendedGrid) (field: Field) (idx: NodeIdx) : single*single =
-    // Fetch surrounding elements
-    let e_n = grid.getElemsSurroundingNode idx
-
-    // Fetch points
-    let p_0 = getVertex grid idx
-    let p_n = e_n |> Array.map (getCentroid grid)
-
-    // Fetch function values
-    let f_n = e_n |> Array.map (Array.get field)
-    let fx_n = f_n |> Array.map fst
-    let fy_n = f_n |> Array.map snd
-
-    // Helper functions
-    let (<->) (a1, a2) (b1, b2) = (a1 - b1, a2 - b2)
-    let inv_norm (px, py) = 1f / sqrt (px*px + py*py)
-
-    // Compute inverse distance weights
-    let w_n = p_n |> Array.map (fun p -> p <-> p_0) |> Array.map inv_norm
-    let W = w_n |> Array.sum
-
-    // Compute mean function value on vertex P_0
-    let fx_0 = (Array.map2 (*) w_n fx_n |> Array.sum) / W
-    let fy_0 = (Array.map2 (*) w_n fy_n |> Array.sum) / W
-    fx_0, fy_0
-
-// Transform centroid (Voronoi) 2D field values to vertex (Delaunay) 2D field values
-let centroidToVertexValues (grid: ExtendedGrid) (field: Field) : Field =
-    let N = ((grid :> IGrid).getVertices ()).Length
-    [| 0..(N-1) |] |> Array.Parallel.map (calcVertexValue grid field)
-
-
-
-
-
-// Cyclic access of array, (A[A.Length] = A[0]), (A[-1] = A[A.Length-1]), etc
-let private access (A: single[]) (idx: int) : single =
-    let l = A.Length
-    let i = idx % l
-    if i >= 0 then A[i] else A[l + i]
-
-// Clough-Tocher interpolation of 1D field inside single triangular cell
-// p_0: position to evaluate, should be within the triangle defined by P_n
-// p_n, f_n, df_n: Positions, values and gradients at cell vertices
-// Ref: C.L. Lawson "C1-Compatible Interpolation Over a Triangle" (1976)
-let private interpolateTriangleCT (p_n: Pos[]) (f_n: single[]) (df_n: Field) (p_0: Pos): single=
-    let idx = [| 0..2 |]
-
-    // Phase 1:
-    let x_n = p_n |> Array.map fst
-    let y_n = p_n |> Array.map snd
-
-    let u = idx |> Array.map (fun i -> (access x_n (i-1)) - (access x_n (i+1)))
-    let v = idx |> Array.map (fun i -> (access y_n (i-1)) - (access y_n (i+1)))
-    let l2 = Array.map2 (fun a b -> a*a + b*b) u v
-
-    let delta = u[0]*v[1] - v[0]*u[1]
-    if (abs delta) < 1.0e-12f then failwith "Co-linear vertices!"
-    let d_inv = 1f / delta
-
-    let x, y = p_0
-    let x_tilde = x - x_n[0]
-    let y_tilde = y - y_n[0]
-
-    let r_1 = d_inv * (u[1]*y_tilde - v[1]*x_tilde)
-    let r_2 = d_inv * (u[2]*y_tilde - v[2]*x_tilde)
-    let r_0 = 1f - (r_1 + r_2)
-    let r = [| r_0; r_1; r_2 |]
-
-    let exterior = r |> Array.exists (fun r_i -> r_i < 0f)
-    if exterior then failwith "Point is outside the triangle!"
-
-    let phi = idx |> Array.map (fun i -> (access r (i+1)) * (access r (i-1)))
-
-    // Phase 2 (Zienkievicz):
-    let rho = idx |> Array.map (fun i ->
-            let r_im1 = access r (i-1)
-            let r_ip1 = access r (i+1)
-            (r[i] * r_ip1*r_ip1 * r_im1*r_im1) / ((1f - r_ip1)*(1f - r_im1))
-        )
-
-    // Phase 3:
-    let fx_n = df_n |> Array.map fst
-    let fy_n = df_n |> Array.map snd
-
-    let h_tilde = idx |> Array.map (fun i ->
-            let fx_ip1 = access fx_n (i+1)
-            let fy_ip1 = access fy_n (i+1)
-            u[i]*fx_ip1 + v[i]*fy_ip1
-        )
-
-    let k_tilde = idx |> Array.map (fun i ->
-            let fx_im1 = access fx_n (i-1)
-            let fy_im1 = access fy_n (i-1)
-            u[i]*fx_im1 + v[i]*fy_im1
-        )
-
-    let g_tilde = idx |> Array.map (fun i ->
-            let r_ip1 = access r (i+1)
-            let r_im1 = access r (i-1)
-            let l2_ip1 = access l2 (i+1)
-            let l2_im1 = access l2 (i-1)
-            let rho_ip1 = access rho (i+1)
-            let rho_im1 = access rho (i-1)
-            (r_ip1 - r_im1) * phi[i] + (3f * rho[i]) * (l2_ip1 - l2_im1) / l2[i] - rho_ip1 + rho_im1
-        )
-
-    let w = idx |> Array.map (fun i ->
-            let f_im1 = access f_n (i-1)
-            let f_ip1 = access f_n (i+1)
-            f_n[i]*r[i] + 0.5f * (h_tilde[i] - k_tilde[i]) * phi[i] +
-            g_tilde[i] * (0.5f * (h_tilde[i] + k_tilde[i]) - f_im1 + f_ip1)
-        )
-    w |> Array.sum
-
-// Clough-Tocher interpolation on vertex 2D field
-let evaluateCloughTocherVertex (grid: ExtendedGrid) (field: Field) (p: Pos) =
-    match grid.tryGetElement p with
-        | None -> None
-        | Some e ->
-            let i0, i1, i2 = getCell grid e
-            let idx = [| i0; i1; i2 |]
-
-            // Split field
-            let U = field |> Array.map fst
-            let V = field |> Array.map snd
-
-            // Position of vertices
-            let p_n = idx |> Array.map (getVertex grid)
-
-            // Field values at vertices
-            let u_n = idx |> Array.map (Array.get U)
-            let v_n = idx |> Array.map (Array.get V)
-
-            // Gradient values at vertices
-            let thrs = 1e-6f
-            let du_n = idx |> Array.map (calcVertexGradient grid U thrs)
-            let dv_n = idx |> Array.map (calcVertexGradient grid V thrs)
-
-            let u1 = interpolateTriangleCT p_n u_n du_n p
-            let v1 = interpolateTriangleCT p_n v_n dv_n p
-            Some (u1, v1)
-
-// Clough-Tocher interpolation on centroid 2D field
-let evaluateCloughTocherCentroid (grid: ExtendedGrid) (field_c: Field) (p: Pos) =
-    let field_v = centroidToVertexValues grid field_c
-    evaluateCloughTocherVertex grid field_v p
-
--- a/src/Fvcom.fs
+++ b/src/Fvcom.fs
@@ -2,206 +2,486 @@ module Oceanbox.FvcomKit.Fvcom

 #nowarn "57"

+open System
 open Microsoft.Research.Science.Data
 open FSharpPlus
-open ProjNet.FSharp
 open Serilog
 open Grid

-type FvcomGrid =
-    {
-        Elem: Elem array
-        Nodes: Node array
-        BBox: BBox
-        Cells: Node array
-        Bathymetry: single []
-        Siglay: single [,]
-        SiglayCenter: single [,]
-        Siglev: single [,]
-        Proj: Projection
-    }
+open Types
+
+type FvcomGrid = {
+    Elem: Elem array
+    Nodes: Node array
+    BBox: BBox
+    Cells: Node array
+    Bathymetry: single[]
+    Siglay: single[,]
+    SiglayCenter: single[,]
+    Siglev: single[,]
+} with
    interface IGrid with
-        member x.getVertex n = x.Nodes[n]
-        member x.getCell n = x.Elem[n]
-        member x.getCellVertices n =
-            let a, b, c = x.Elem[n]
-            x.Nodes[a], x.Nodes[b], x.Nodes[c]
-        member x.getVertices() = x.Nodes
-        member x.getCells() = x.Elem
-        member x.getBoundingBox() = x.BBox
-        member x.Projection = x.Proj
-    static member empty =
-        {
-            Elem = Array.empty
-            Nodes = Array.empty
-            BBox = BBox.empty
-            Cells = Array.empty
-            Bathymetry = Array.empty
-            Siglay = Array2D.zeroCreate 0 0
-            SiglayCenter = Array2D.zeroCreate 0 0
-            Siglev = Array2D.zeroCreate 0 0
-            Proj = WebMercator
-        }
-    member this.ToGrid() =
-        {
-            Elem = this.Elem
-            Nodes = this.Nodes
-            BBox = this.BBox
-            Proj = this.Proj
-        }
+        member this.getVertex n = this.Nodes[n]
+        member this.getCell n = this.Elem[n]
+        member this.getCellVertices n =
+            let a, b, c = this.Elem[n]
+            this.Nodes[a], this.Nodes[b], this.Nodes[c]
+        member this.getVertices() = this.Nodes
+        member this.getCells() = this.Elem
+        member this.getBoundingBox() = this.BBox
+
+    static member empty = {
+        Elem = Array.empty
+        Nodes = Array.empty
+        BBox = BBox.empty
+        Cells = Array.empty
+        Bathymetry = Array.empty
+        Siglay = Array2D.zeroCreate 0 0
+        SiglayCenter = Array2D.zeroCreate 0 0
+        Siglev = Array2D.zeroCreate 0 0
+    }
+
+    member this.ToGrid() = { Elem = this.Elem; Nodes = this.Nodes; BBox = this.BBox }

 let getNumFrames (ds: DataSet) =
    try
        ds.Dimensions["time"].Length
-    with
-    | e ->
+    with e ->
        Log.Error $"{e}"
        0

+let getNumSiglay (ds: DataSet) =
+    try
+        ds.Dimensions["siglay"].Length
+    with e ->
+        Log.Error $"{e}"
+        0
+
+let getNumSiglev (ds: DataSet) =
+    try
+        ds.Dimensions["siglev"].Length
+    with e ->
+        Log.Error $"{e}"
+        0
+
+let getTime (ds: DataSet) n =
+    try
+        let ts = ds["Times"].GetData () :?> byte[,]
+
+        ts[n, *] |> Array.map char |> System.String |> System.DateTime.Parse |> Some
+    with e ->
+        Log.Error $"getTime exception: {e.Message}"
+        None
+
+let getTimeSpanSinceStart (ds: DataSet) n =
+    try
+        let days = ds["Itime"].GetData () :?> int[]
+        let msec = ds["Itime2"].GetData () :?> int[]
+        let t0 = TimeSpan.FromDays days[n]
+        let t1 = TimeSpan.FromMilliseconds (float msec[n])
+        t0 + t1 |> Some
+    with e ->
+        Log.Error $"getTimeInDays exception: {e.Message}"
+        None
+
+let readTauc (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let tauc = ds["tauc"].GetData ([| t; 0 |], [| 1; n |]) :?> single[,]
+        tauc[0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
 let readUV (ds: DataSet) t l =
    try
        let n = ds.Dimensions["nele"].Length
-        let u = ds[ "u" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
-        let v = ds[ "v" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
+        let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
        Array.zip u[0, 0, *] v[0, 0, *]
-    with
-    | e ->
+    with e ->
        Log.Error $"{e}"
        Array.empty

 let readUV' (ds: DataSet) t l =
-    readUV ds t l
-    |> Array.collect (fun (x, y) -> [| x; y |])
+    readUV ds t l |> Array.collect (fun (x, y) -> [| x; y |])
+
+let readUVs (ds: DataSet) t l es =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let u' = es |> Array.map (fun i -> u[0, 0, i])
+        let v' = es |> Array.map (fun i -> v[0, 0, i])
+        Array.zip u' v'
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+/// <summary>
+/// Reads properties 'u' and 'v' from FVCOM dataset at the given time and depth.
+/// </summary>
+/// <param name="ds">The NetCDF dataset to open</param>
+/// <param name="t">The given time frame</param>
+/// <param name="l">The layer of depth from where to read the values</param>
+/// <param name="e0">The starting index of the uv range</param>
+/// <param name="en">The ending index</param>
+/// <returns>Array of tuples, and empty on any errors.</returns>
+let readUVRange (ds: DataSet) t l e0 en =
+    try
+        let u = ds["u"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
+        let v = ds["v"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
+        Array.zip u[0, 0, *] v[0, 0, *]
+    with e ->
+        Log.Error (
+            e,
+            "FvcomKit.Fvcom.readUVRange exception with input: time {Time}, layer {Layer}, starting elem {Elem0} and ending elem {ElemN}",
+            t,
+            l,
+            e0,
+            en
+        )
+        Array.empty
+
+let readOmega (ds: DataSet) t l =
+    try
+        let n = ds.Dimensions["node"].Length
+        let omega = ds["omega"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        omega[0, 0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readOmegaBlock (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["node"].Length
+        let l = ds.Dimensions["siglev"].Length
+        let u = ds["omega"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
+        u[0, *, *]
+    with e ->
+        Log.Error $"{e}"
+        Array2D.zeroCreate 0 0
+
+let readOmegas (ds: DataSet) t l es =
+    try
+        let n = ds.Dimensions["node"].Length
+        let omega = ds["omega"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        es |> Array.map (fun i -> omega[0, 0, i])
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readU (ds: DataSet) t l =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        u[0, 0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readUBlock (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let l = ds.Dimensions["siglay"].Length
+        let u = ds["u"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
+        u[0, *, *]
+    with e ->
+        Log.Error $"{e}"
+        Array2D.zeroCreate 0 0
+
+let readV (ds: DataSet) t l =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        v[0, 0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readVBlock (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let l = ds.Dimensions["siglay"].Length
+        let v = ds["v"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
+        v[0, *, *]
+    with e ->
+        Log.Error $"{e}"
+        Array2D.zeroCreate 0 0
+
+let readWw (ds: DataSet) t l =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let ww = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        ww[0, 0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readWwBlock (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let l = ds.Dimensions["siglay"].Length
+        let w = ds["ww"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
+        w[0, *, *]
+    with e ->
+        Log.Error $"{e}"
+        Array2D.zeroCreate 0 0
+
+let readWws (ds: DataSet) t l es =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let ww = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        es |> Array.map (fun i -> ww[0, 0, i])
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readUVW (ds: DataSet) t l =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let w = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        Array.zip3 u[0, 0, *] v[0, 0, *] w[0, 0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readUVWs (ds: DataSet) t l es =
+    try
+        let n = ds.Dimensions["nele"].Length
+        let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let w = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
+        let u' = es |> Array.map (fun i -> u[0, 0, i])
+        let v' = es |> Array.map (fun i -> v[0, 0, i])
+        let w' = es |> Array.map (fun i -> w[0, 0, i])
+        Array.zip3 u' v' w'
+    with e ->
+        Log.Error $"{e}"
+        Array.empty
+
+let readUVWRange (ds: DataSet) t l e0 en =
+    try
+        let u = ds["u"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
+        let v = ds["v"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
+        let w = ds["ww"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
+        Array.zip3 u[0, 0, *] v[0, 0, *] w[0, 0, *]
+    with e ->
+        Log.Error $"{e}"
+        Array.empty

 let readTemp (ds: DataSet) t l =
    try
        let n = ds.Dimensions["node"].Length
-        let temp =
-            ds[ "temp" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
+        let temp = ds["temp"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
        temp[0, 0, *]
-    with
-    | e ->
+    with e ->
        Log.Error $"{e}"
        Array.empty

+let readTempBlock (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["node"].Length
+        let l = ds.Dimensions["siglay"].Length
+        let sal = ds["temp"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
+        sal[0, *, *]
+    with e ->
+        Log.Error $"{e}"
+        Array2D.zeroCreate 0 0
+
 let readSalinity (ds: DataSet) t l =
    try
        let n = ds.Dimensions["node"].Length
-        let sal =
-            ds["salinity"]
-                .GetData([| t; l; 0 |], [| 1; 1; n |])
-            :?> single [,,]
+        let sal = ds["salinity"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
        sal[0, 0, *]
-    with
-    | e ->
+    with e ->
        Log.Error $"{e}"
        Array.empty

+let readSalinityBlock (ds: DataSet) t =
+    try
+        let n = ds.Dimensions["node"].Length
+        let l = ds.Dimensions["siglay"].Length
+        let sal = ds["salinity"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
+        sal[0, *, *]
+    with e ->
+        Log.Error $"{e}"
+        Array2D.zeroCreate 0 0
+
+
 let readArt1 (ds: DataSet) =
    try
-        ds[ "art1" ].GetData() :?> single []
-    with
-    | e ->
+        ds["art1"].GetData () :?> single[]
+    with e ->
        Log.Error $"{e}"
        Array.empty

 let readZeta (ds: DataSet) t =
    try
        let n = ds.Dimensions["node"].Length
-        let zeta = ds[ "zeta" ].GetData([| t; 0 |], [| 1; n |]) :?> single [,]
+        let zeta = ds["zeta"].GetData ([| t; 0 |], [| 1; n |]) :?> single[,]
        zeta[0, *]
-    with
-    | e ->
+    with e ->
        Log.Error $"{e}"
        Array.empty

 let readBathymetry (ds: DataSet) =
    try
-        ds[ "h" ].GetData() :?> single []
-    with
-    | e ->
+        ds["h"].GetData () :?> single[]
+    with e ->
+        Log.Error e.Message
+        Array.empty
+
+let tryReadBathymetry (ds: DataSet) =
+    if ds.Variables.Contains "h" then
+        ds["h"].GetData () :?> single[] |> Some
+    else
+        None
+
+let readBathymetryAtCenters (ds: DataSet) =
+    let h = readBathymetry ds
+    let lc = ds.Dimensions["nele"].Length
+    try
+        let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
+        [| 0 .. lc - 1 |]
+        |> Array.map (fun i -> (h[nv[0, i]] + h[nv[1, i]] + h[nv[2, i]]) / 3.0f)
+    with e ->
        Log.Error e.Message
        Array.empty

 let readBathymetryAtCenter (ds: DataSet) e =
    try
-        let h = ds[ "h_center" ].GetData() :?> single []
+        let h = readBathymetryAtCenters ds
        h[e]
-    with
-    | e ->
+    with e ->
        Log.Error e.Message
        -1f

 let readSiglev (ds: DataSet) n =
    try
-        let siglev = ds[ "siglev" ].GetData() :?> single [,]
-        siglev[*, n]
-    with
-    | err ->
+        let l = ds.Dimensions["siglev"].Length
+        let siglev = ds["siglev"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
+        siglev[*, 0]
+    with err ->
        Log.Error $"{err}"
        Array.empty

+let tryReadSiglev (ds: DataSet) n =
+    if ds.Variables.Contains "siglev" then
+        let l = ds.Dimensions["siglev"].Length
+        let siglev = ds["siglev"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
+        Some siglev[*, 0]
+    else
+        None
+
 let readSiglay (ds: DataSet) n =
    try
-        let siglay = ds[ "siglay" ].GetData() :?> single [,]
-        siglay[*, n]
-    with
-    | err ->
+        let l = ds.Dimensions["siglay"].Length
+        let siglay = ds["siglay"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
+        siglay[*, 0]
+    with err ->
        Log.Error $"{err}"
        Array.empty

-let readSiglayCenter (ds: DataSet) n =
-    try
-        let siglay = ds[ "siglay_center" ].GetData() :?> single [,]
-        siglay[*, n]
-    with
-    | err ->
-        Log.Error $"{err}"
-        Array.empty
+let tryReadSiglay (ds: DataSet) n =
+    if ds.Variables.Contains "siglay" then
+        let l = ds.Dimensions["siglay"].Length
+        let siglay = ds["siglay"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
+        Some siglay[*, 0]
+    else
+        None

 let readSiglevAtCenter (ds: DataSet) e =
    try
-        let siglev = ds[ "siglev_center" ].GetData() :?> single [,]
-        siglev[*, e]
-    with
-    | e ->
+        let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
+        let s1 = readSiglev ds nv[0, e]
+        let s2 = readSiglev ds nv[1, e]
+        let s3 = readSiglev ds nv[2, e]
+        [| 0 .. s1.Length - 1 |] |> Array.map (fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
+    with e ->
        Log.Error $"{e}"
        Array.empty

 let readSiglayAtCenter (ds: DataSet) e =
    try
-        let siglay = ds[ "siglay_center" ].GetData() :?> single [,]
-        siglay[*, e]
-    with
-    | e ->
+        let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
+        let s1 = readSiglay ds nv[0, e]
+        let s2 = readSiglay ds nv[1, e]
+        let s3 = readSiglay ds nv[2, e]
+        [| 0 .. s1.Length - 1 |] |> Array.map (fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
+    with e ->
        Log.Error $"{e}"
        Array.empty

+let tryReadSiglayAtCenter (ds: DataSet) e =
+    if ds.Variables.Contains "siglay" then
+        let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
+        let s1 = readSiglay ds nv[0, e]
+        let s2 = readSiglay ds nv[1, e]
+        let s3 = readSiglay ds nv[2, e]
+        [| 0 .. s1.Length - 1 |]
+        |> Array.map (fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
+        |> Some
+    else
+        None
+
+// deprecated
+let readSiglayCenter = readSiglayAtCenter
+
 module Siglay =
    let readSiglay (ds: DataSet) =
        try
-            ds[ "siglay" ].GetData() :?> single [,]
-        with
-        | err ->
+            ds["siglay"].GetData () :?> single[,]
+        with err ->
            Log.Error $"{err}"
            Array2D.zeroCreate 0 0

+    let tryReadSiglay (ds: DataSet) =
+        if ds.Variables.Contains "siglay" then
+            ds["siglay"].GetData () :?> single[,] |> Some
+        else
+            None
+
    let readSiglayAtCenter (ds: DataSet) =
+        let siglay = readSiglay ds
        try
-            ds[ "siglay_center" ].GetData() :?> single [,]
-        with
-        | e ->
+            let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
+            let l1 = Array2D.length1 siglay
+            let l2 = Array2D.length2 nv
+            let sc = Array2D.zeroCreate l1 l2
+            sc
+            |> Array2D.mapi (fun i j _ -> (siglay[i, nv[0, j]] + siglay[i, nv[1, j]] + siglay[i, nv[2, j]]) / 3.0f)
+        with e ->
            Log.Error $"{e}"
            Array2D.zeroCreate 0 0

+    let tryReadSiglayAtCenter (ds: DataSet) =
+        if ds.Variables.Contains "siglay" then
+            readSiglayAtCenter ds |> Some
+        else
+            None
+
    let readUv (ds: DataSet) e t =
        try
            let l = ds.Dimensions["siglay"].Length
-            let u = ds[ "u" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
-            let v = ds[ "v" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
+            let u = ds["u"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
+            let v = ds["v"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
            Array.zip u[0, *, 0] v[0, *, 0]
-        with
-        | err ->
+        with err ->
+            Log.Warning $"readUv {e} {t}"
+            Log.Error $"{err}"
+            Array.empty
+
+    let readUvw (ds: DataSet) e t =
+        try
+            let l = ds.Dimensions["siglay"].Length
+            let u = ds["u"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
+            let v = ds["v"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
+            let w = ds["ww"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
+            Array.zip3 u[0, *, 0] v[0, *, 0] w[0, *, 0]
+        with err ->
            Log.Warning $"readUv {e} {t}"
            Log.Error $"{err}"
            Array.empty
@@ -209,12 +489,10 @@ module Siglay =
    let readUv' (ds: DataSet) e t =
        try
            let l = ds.Dimensions["siglay"].Length
-            let u = ds[ "u" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
-            let v = ds[ "v" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
-            Array.zip u[0, *, 0] v[0, *, 0]
-            |> Array.collect (fun (x, y) -> [| x; y |])
-        with
-        | err ->
+            let u = ds["u"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
+            let v = ds["v"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
+            Array.zip u[0, *, 0] v[0, *, 0] |> Array.collect (fun (x, y) -> [| x; y |])
+        with err ->
            Log.Warning $"readUv' {e} {t}"
            Log.Error $"{err}"
            Array.empty
@@ -222,10 +500,9 @@ module Siglay =
    let readTemp (ds: DataSet) n t =
        try
            let l = ds.Dimensions["siglay"].Length
-            let t = ds[ "temp" ].GetData([| t; 0; n |], [| 1; l; 1 |]) :?> single [,,]
+            let t = ds["temp"].GetData ([| t; 0; n |], [| 1; l; 1 |]) :?> single[,,]
            t[0, *, 0]
-        with
-        | err ->
+        with err ->
            Log.Warning $"readTemp {n} {t}"
            Log.Error $"{err}"
            Array.empty
@@ -233,60 +510,100 @@ module Siglay =
    let readSalinity (ds: DataSet) n t =
        try
            let l = ds.Dimensions["siglay"].Length
-            let t =
-                ds["salinity"]
-                    .GetData([| t; 0; n |], [| 1; l; 1 |])
-                :?> single [,,]
+            let t = ds["salinity"].GetData ([| t; 0; n |], [| 1; l; 1 |]) :?> single[,,]
            t[0, *, 0]
-        with
-        | err ->
+        with err ->
            Log.Warning $"readS {n} {t}"
            Log.Error $"{err}"
            Array.empty

+module Singular =
+    let readUv (ds: DataSet) e t d =
+        try
+            let u = ds["u"].GetData ([| t; d; e |], [| 1; 1; 1 |]) :?> single[,,]
+            let v = ds["v"].GetData ([| t; d; e |], [| 1; 1; 1 |]) :?> single[,,]
+            u[0, 0, 0], v[0, 0, 0]
+        with err ->
+            Log.Warning $"readUv {e} {t}"
+            Log.Error $"{err}"
+            0f, 0f
+
+    let readTemp (ds: DataSet) n t d : single =
+        try
+            let t = ds["temp"].GetData ([| t; d; n |], [| 1; 1; 1 |]) :?> single[,,]
+            t[0, 0, 0]
+        with err ->
+            Log.Warning $"readTemp {n} {t}"
+            Log.Error $"{err}"
+            0f
+
+    let readSalinity (ds: DataSet) n t d =
+        try
+            let t = ds["salinity"].GetData ([| t; d; n |], [| 1; 1; 1 |]) :?> single[,,]
+            t[0, 0, 0]
+        with err ->
+            Log.Warning $"readS {n} {t}"
+            Log.Error $"{err}"
+            0f
+
+    let readZeta (ds: DataSet) n t =
+        try
+            let zeta = ds["zeta"].GetData ([| t; n |], [| 1; 1 |]) :?> single[,]
+            zeta[0, 0]
+        with e ->
+            Log.Error $"{e}"
+            0f
+
 let getBBox (ds: DataSet) : BBox =
    try
-        let x = ds[ "x" ].GetData() :?> single []
-        let y = ds[ "y" ].GetData() :?> single []
+        let x = ds["x"].GetData () :?> single[]
+        let y = ds["y"].GetData () :?> single[]
        let minX = Array.min x
        let maxX = Array.max x
        let minY = Array.min y
        let maxY = Array.max y
-        let center = minX + (maxX - minX) / 2f, minY + (maxY - minY) / 2f
+        let center = float (minX + maxX) / 2., float (minY + maxY) / 2.
        {
-            minX = minX
-            maxX = maxX
-            minY = minY
-            maxY = maxY
+            minX = float minX
+            maxX = float maxX
+            minY = float minY
+            maxY = float maxY
            center = center
        }
-    with
-    | e ->
+    with e ->
        Log.Error $"{e}"
        BBox.empty

 let getGrid (ds: DataSet) : FvcomGrid =
    try
-        let x = ds[ "x" ].GetData() :?> single []
-        let y = ds[ "y" ].GetData() :?> single []
-        let xc = ds[ "xc" ].GetData() :?> single []
-        let yc = ds[ "yc" ].GetData() :?> single []
-        let nv = ds[ "nv" ].GetData() :?> int [,]
+        let x = ds["x"].GetData () :?> single[] |> Array.map float
+        let y = ds["y"].GetData () :?> single[] |> Array.map float
+        let xc = ds["xc"].GetData () :?> single[] |> Array.map float
+        let yc = ds["yc"].GetData () :?> single[] |> Array.map float
+        let nv = ds["nv"].GetData () :?> int[,]

-        let h = readBathymetry ds
-        let siglay = Siglay.readSiglay ds
-        let siglay_c = Siglay.readSiglayAtCenter ds
-        let siglev = ds[ "siglev" ].GetData() :?> single [,]
+        let h =
+            match tryReadBathymetry ds with
+            | Some b -> b
+            | None -> Array.empty
+        let siglay =
+            match Siglay.tryReadSiglay ds with
+            | Some s -> s
+            | None -> Array2D.zeroCreate 0 0
+        let siglay_c =
+            match Siglay.tryReadSiglayAtCenter ds with
+            | Some s -> s
+            | None -> Array2D.zeroCreate 0 0
+        let siglev =
+            if ds.Variables.Contains "siglev" then
+                ds["siglev"].GetData () :?> single[,]
+            else
+                Array2D.zeroCreate 0 0
        let elem =
            Array.zip3 nv[0, *] nv[1, *] nv[2, *]
            |> Array.map (fun (a, b, c) -> a - 1, b - 1, c - 1)
        let nds = Array.zip x y
        let cells = Array.zip xc yc
-        let proj =
-            try
-                ds.Metadata["CoordinateProjection"] :?> string
-            with
-            | _ -> ""
        {
            Elem = elem
            Nodes = nds
@@ -296,26 +613,65 @@ let getGrid (ds: DataSet) : FvcomGrid =
            Siglay = siglay
            SiglayCenter = siglay_c
            Siglev = siglev
-            Proj = Projection.FromString proj
        }
-    with
-    | e ->
+    with e ->
+        Log.Error $"{e}"
+        FvcomGrid.empty
+
+let getGridLonLat (ds: DataSet) : FvcomGrid =
+    try
+        let x = ds["lon"].GetData () :?> single[] |> Array.map float
+        let y = ds["lat"].GetData () :?> single[] |> Array.map float
+        let xc = ds["lonc"].GetData () :?> single[] |> Array.map float
+        let yc = ds["latc"].GetData () :?> single[] |> Array.map float
+        let nv = ds["nv"].GetData () :?> int[,]
+
+        let h =
+            match tryReadBathymetry ds with
+            | Some b -> b
+            | None -> Array.empty
+        let siglay =
+            match Siglay.tryReadSiglay ds with
+            | Some s -> s
+            | None -> Array2D.zeroCreate 0 0
+        let siglay_c =
+            match Siglay.tryReadSiglayAtCenter ds with
+            | Some s -> s
+            | None -> Array2D.zeroCreate 0 0
+        let siglev =
+            if ds.Variables.Contains "siglev" then
+                ds["siglev"].GetData () :?> single[,]
+            else
+                Array2D.zeroCreate 0 0
+        let elem =
+            Array.zip3 nv[0, *] nv[1, *] nv[2, *]
+            |> Array.map (fun (a, b, c) -> a - 1, b - 1, c - 1)
+        let nds = Array.zip x y
+        let cells = Array.zip xc yc
+        {
+            Elem = elem
+            Nodes = nds
+            BBox = getBBox ds
+            Cells = cells
+            Bathymetry = h
+            Siglay = siglay
+            SiglayCenter = siglay_c
+            Siglev = siglev
+        }
+    with e ->
        Log.Error $"{e}"
        FvcomGrid.empty

 let getNbve (ds: DataSet) =
-    let nbve = ds[ "nbve" ].GetData() :?> int [,]
+    let nbve = ds["nbve"].GetData () :?> int[,]
    [|
        for i = 0 to (Array2D.length2 nbve - 1) do
-            nbve[*, i]
-            |> Array.filter ((<>) 0)
-            |> Array.map (fun x -> x - 1)
-            |> Array.rev
+            nbve[*, i] |> Array.filter ((<>) 0) |> Array.map (fun x -> x - 1) |> Array.rev
    |]

-let projectFvcomGrid proj (grid: FvcomGrid) : FvcomGrid =
-    { grid with
+let projectFvcomGrid proj (grid: FvcomGrid) : FvcomGrid = {
+    grid with
        Nodes = grid.Nodes |> Array.Parallel.map proj
        BBox = projectBBox proj grid.BBox
        Cells = grid.Cells |> Array.Parallel.map proj
-    }
+}
--- a/src/Gradient.fs
+++ b/src/Gradient.fs
@@ -1,108 +0,0 @@
-module Oceanbox.FvcomKit.Gradient
-
-open Grid
-
-let private getCell (grid: ExtendedGrid) (idx: ElemIdx) =
-    (grid :> IGrid).getCell idx
-
-let private getVertex (grid: ExtendedGrid) (idx: NodeIdx) =
-    (grid :> IGrid).getVertex idx
-
-// Compute contributions to Q and s from a single neighboring point, given current estimate of gradient df
-let private calcContributionQandS (df: single*single) (p0: Pos) (f0: single) (p1: Pos) (f1: single) =
-    let x0, y0 = p0
-    let x1, y1 = p1
-
-    let ex = x1 - x0
-    let ey = y1 - y0
-
-    let L = sqrt (ex*ex + ey*ey)
-    let L3 = L*L*L
-
-    let dfx, dfy = df
-    let df2 = -ex*dfx - ey*dfy
-
-    let q1 = 4f * ex * ex / L3
-    let q2 = 4f * ex * ey / L3
-    let q3 = 4f * ey * ey / L3
-
-    let s1 = (6f*(f0 - f1) - 2f*df2) * ex / L3
-    let s2 = (6f*(f0 - f1) - 2f*df2) * ey / L3
-
-    let Q = [| q1; q2; q2; q3 |]
-    let s = [| s1; s2 |]
-    Q, s
-
-// Set up linear system to solve (Qx = s), accumulated from N neighboring points
-let private setupQandS (p0: Pos) (f0: single) (p_n: Pos[]) (f_n: single[]) (df: single*single) =
-    let calcQS_n = calcContributionQandS df p0 f0
-
-    let (<+>) a b = Array.map2 (+) a b
-    let rec accumulateQandS Q s idx =
-        match idx with
-        | [] -> Q, s
-        | x::xs ->
-            let Q', s' = calcQS_n p_n[x] f_n[x]
-            accumulateQandS (Q <+> Q') (s <+> s') xs
-
-    let Q0 = [| 0f; 0f; 0f; 0f |]
-    let s0 = [| 0f; 0f |]
-    let idx = [ 0..(p_n.Length-1) ]
-    accumulateQandS Q0 s0 idx
-
-// Estimates gradient in point P_0 based on neighboring function values
-// p0: Point of interest
-// f0: Function value in POI
-// p_n: Neighboring points
-// f_n: Function values in neighboring points
-// thrs: Convergence threshold in iterative procedure
-//
-// Ref: Based on SciPy implementation which uses the following
-// G. Nielson "A method for interpolating scattered data based upon a minimum norm network" (1983)
-// R. J. Renka and A. K. Cline "A Triangle-based C1 interpolation method" (1984)
-let calcGradientFromNeighbors (p0: Pos) (f0: single) (p_n: Pos[]) (f_n: single[]) (thrs: single) : single*single =
-    let setupQS = setupQandS p0 f0 p_n f_n
-
-    let rec convergeGradient (error: single) (df: single*single) =
-        if error < thrs then
-            df
-        else
-            let Q, s = setupQS df
-            let detQ = Q[0]*Q[3] - Q[1]*Q[2]
-
-            if (abs detQ) < thrs*thrs then
-                df
-            else
-                let r0 = ( Q[3]*s[0] - Q[1]*s[1]) / detQ
-                let r1 = (-Q[2]*s[0] + Q[0]*s[1]) / detQ
-                let dfx, dfy = df
-                let change = max (abs (dfx + r0)) (abs (dfy + r1))
-                convergeGradient change (-r0, -r1)
-
-    let err0 = 1f
-    let grad0 = 0f, 0f
-    convergeGradient err0 grad0
-
-// Given vertex 1D field values, compute gradient in single vertex
-let calcVertexGradient (grid: ExtendedGrid) (field: single[]) (thrs: single) (idx: NodeIdx) : single*single =
-    let neighbors =
-        grid.getNodesSurroundingNode idx
-        |> Array.filter (fun i -> i <> idx)
-
-    // Get data from point of interest
-    let p0 = getVertex grid idx
-    let f0 = field[idx]
-
-    // Get data from neighbors of POI
-    let p_n = neighbors |> Array.map (getVertex grid)
-    let f_n = neighbors |> Array.map (Array.get field)
-
-    calcGradientFromNeighbors p0 f0 p_n f_n thrs
-
-// Given vertex 1D field values, compute gradients on entire grid
-let calcVertexGradients (grid: ExtendedGrid) (field: single[]) (thrs: single) : Field =
-    let N = ((grid :> IGrid).getVertices ()).Length
-    [| 0..(N-1) |] |> Array.map (calcVertexGradient grid field thrs)
-
-
-
--- a/src/Grid.fs
+++ b/src/Grid.fs
@@ -6,117 +6,90 @@ open FSharpPlus.Control
 open ProjNet.CoordinateSystems
 open Serilog
 open ProjNet.FSharp
+open MessagePack
 open MBrace.FsPickler
-open FsKDTree
-open KdTree // C# version
+//open FsKDTree
+open KdTree // NOTE: C# version
+
+open Types

 type NodeIdx = int
 type ElemIdx = int
 type Edge = int * int

 type Elem = NodeIdx * NodeIdx * NodeIdx
-type Node = single * single
+type Node = float * float

-type Pos = single * single
-type Field = (single*single) []
-
-type BBox =
-    {
-        minX: single
-        maxX: single
-        minY: single
-        maxY: single
-        center: single * single
-    }
-    static member empty =
-        {
-            minX = Single.MaxValue
-            maxX = Single.MinValue
-            minY = Single.MaxValue
-            maxY = Single.MinValue
-            center = 0f, 0f
-        }
+type Pos = float * float
+type Leaf<'a> = { Pos: Pos; Data: 'a }
+type Field = (float * float) array

 type Cell = NodeIdx * NodeIdx * NodeIdx
-type Vertex = single * single

 type IGrid =
    abstract getVertex: int -> Vertex
    abstract getCell: int -> Cell
    abstract getCellVertices: int -> Vertex * Vertex * Vertex
-    abstract getVertices: unit -> Vertex []
-    abstract getCells: unit -> Cell []
+    abstract getVertices: unit -> Vertex array
+    abstract getCells: unit -> Cell array
    abstract getBoundingBox: unit -> BBox
-    abstract Projection: Projection

-type Grid =
-    {
-        Elem: Elem array
-        Nodes: Node array
-        BBox: BBox
-        Proj: Projection
-    }
+type Grid = {
+    Elem: Elem array
+    Nodes: Node array
+    BBox: BBox
+} with
    interface IGrid with
-        member x.getVertex n = x.Nodes[n]
-        member x.getCell n = x.Elem[n]
-        member x.getCellVertices n =
-            let a, b, c = x.Elem[n]
-            x.Nodes[a], x.Nodes[b], x.Nodes[c]
-        member x.getVertices() = x.Nodes
-        member x.getCells() = x.Elem
-        member x.getBoundingBox() = x.BBox
-        member x.Projection = x.Proj
-    static member empty =
-        {
-            Elem = Array.empty
-            Nodes = Array.empty
-            BBox = BBox.empty
-            Proj = WebMercator
-        }
+        member this.getVertex n = this.Nodes[n]
+        member this.getCell n = this.Elem[n]
+        member this.getCellVertices n =
+            let a, b, c = this.Elem[n]
+            this.Nodes[a], this.Nodes[b], this.Nodes[c]
+        member this.getVertices() = this.Nodes
+        member this.getCells() = this.Elem
+        member this.getBoundingBox() = this.BBox
+    static member empty = { Elem = Array.empty; Nodes = Array.empty; BBox = BBox.empty }

-type ElemsAroundNode = Map<NodeIdx, ElemIdx []>
-type NodesAroundNode = Map<NodeIdx, NodeIdx []>
+type ElemsAroundNode = Map<NodeIdx, ElemIdx array>
+type NodesAroundNode = Map<NodeIdx, NodeIdx array>
+type ElemsAroundElem = Map<ElemIdx, ElemIdx array>

-type NeighborIndex =
-    {
-        ElemsAroundNode: ElemsAroundNode
-        NodesAroundNode: NodesAroundNode
-    }
-    static member empty = { ElemsAroundNode = Map.empty; NodesAroundNode = Map.empty }
-
-type private Ean = Map<NodeIdx, ElemIdx list>
+type NeighborIndex = {
+    ElemsAroundNode: ElemsAroundNode
+    NodesAroundNode: NodesAroundNode
+    ElemsAroundElem: ElemsAroundElem
+} with
+    static member empty = { ElemsAroundNode = Map.empty; NodesAroundNode = Map.empty; ElemsAroundElem = Map.empty }

 // NOTE(SimenLK): The amount of items to be stored in the trees leafs
-let treeLeafSize = LeafNodeSize 64
+// let treeLeafSize = LeafNodeSize 64

-let private createTree (points: Leaf<single, int>[]) =
-    let tree = KdTree<single, int>(2, KdTree.Math.FloatMath())
-    points
-    |> Array.iter (fun a -> tree.Add([| a.Pos.X; a.Pos.Y |], a.Data) |> ignore)
+let private createTree (points: Leaf<int> array) =
+    let tree = KdTree<float, int> (2, KdTree.Math.DoubleMath ())
+
+    do points |> Array.iter (fun a -> tree.Add ([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
    if points.Length > 0 then
-        tree.Balance()
+        do tree.Balance ()
    else
-        Log.Warning $"Empty kd-tree"
+        do Log.Warning $"Empty kd-tree"
+
    tree

+type private Ean = Map<NodeIdx, ElemIdx array>

 let private makeElemsSurroundingNodeMap (elem: Elem array) : ElemsAroundNode =
-    let addElIdx k v (nodes: Ean) =
-        Map.tryFind k nodes
-        |> Option.defaultWith (fun () -> [])
-        |> fun nds -> Map.add k (v :: nds) nodes
-    elem
-    |> Array.fold
-        (fun (n, acc) (a, b, c) ->
-            let acc' =
-                acc
-                |> addElIdx a n
-                |> addElIdx b n
-                |> addElIdx c n
-            n + 1, acc')
-        (0, Map.empty)
+    let addElIdx k v (nodes: Ean) : Ean =
+        let nds =
+            Map.tryFind k nodes
+            |> Option.defaultValue [||]
+        nodes |> Map.add k (Array.append [|v|] nds)
+    let folder (n, acc) (a, b, c) =
+        let acc' = acc |> addElIdx a n |> addElIdx b n |> addElIdx c n
+        n + 1, acc'
+
+    ((0, Map.empty), elem)
+    ||> Array.fold folder
    |> snd
-    |> Map.mapValues toArray

 let private makeElemsSurroundingNodeMap' (elem: Elem array) : ElemsAroundNode =
    let addElemIdx k v nodes =
@@ -127,12 +100,9 @@ let private makeElemsSurroundingNodeMap' (elem: Elem array) : ElemsAroundNode =
    elem
    |> Array.fold
        (fun (n, acc) (a, b, c) ->
-            let acc' =
-                acc
-                |> addElemIdx a n
-                |> addElemIdx b n
-                |> addElemIdx c n
-            n + 1, acc')
+            let acc' = acc |> addElemIdx a n |> addElemIdx b n |> addElemIdx c n
+            n + 1, acc'
+        )
        (0, Map.empty)
    |> snd
    |> Map.mapValues toArray
@@ -143,34 +113,51 @@ let private makeNodesSurroudingNodeMap (n2e: ElemsAroundNode) (elem: Elem array)
        n
        |> Array.collect (fun x ->
            let n1, n2, n3 = elem[x]
-            [| n1; n2; n3 |])
-        |> Array.distinct)
+            [| n1; n2; n3 |]
+        )
+        |> Array.distinct
+    )

-let getSurrounding (idx: Map<int, int []>) (a, b, c) =
-    [| idx[a]; idx[b]; idx[c] |]
-    |> Array.concat
-    |> Array.distinct
+let private makeElemsSurroundingElemMap (ean: ElemsAroundNode) (elem: Elem array) : ElemsAroundElem =
+    elem
+    |> Array.mapi (fun elemIdx (n1, n2, n3) ->
+        // For each element, find all elements that share any of its nodes
+        let surroundingElems =
+            [| ean[n1]; ean[n2]; ean[n3] |]
+            |> Array.concat
+            |> Array.distinct
+            |> Array.filter (fun x -> x <> elemIdx) // Remove self
+        elemIdx, surroundingElems
+    )
+    |> Map.ofArray
+
+let getSurrounding (idx: Map<int, int[]>) (a, b, c) =
+    [| idx[a]; idx[b]; idx[c] |] |> Array.concat |> Array.distinct

 let makeNeighborIndex (grid: IGrid) =
    let elem = grid.getCells ()
    let ean = makeElemsSurroundingNodeMap elem
+
    {
        ElemsAroundNode = ean
        NodesAroundNode = makeNodesSurroudingNodeMap ean elem
+        ElemsAroundElem = makeElemsSurroundingElemMap ean elem
    }

 let getElemsSurroundingNode (idx: NeighborIndex) n = idx.ElemsAroundNode[n]

 let getNodesSurroundingNode (idx: NeighborIndex) n = idx.NodesAroundNode[n]

-let getNodesSurroundingElem (idx: NeighborIndex) (grid: Grid) e =
-    getSurrounding idx.NodesAroundNode grid.Elem[e]
+let getNodesSurroundingElem (idx: NeighborIndex) e = idx.ElemsAroundElem[e]
+
+// let getNodesSurroundingElem (idx: NeighborIndex) (grid: Grid) e =
+    // getSurrounding idx.NodesAroundNode grid.Elem[e]

 let getElemsSurroundingElem (idx: NeighborIndex) (grid: Grid) e =
    getSurrounding idx.ElemsAroundNode grid.Elem[e]

 let calcBBoxCenter (x: BBox) =
-    x.minX + (x.maxX - x.minX) / 2.0f, x.minY + (x.maxY - x.minY) / 2.0f
+    x.minX + (x.maxX - x.minX) / 2.0, x.minY + (x.maxY - x.minY) / 2.0

 let calcBBox nodes =
    let x, y = Array.unzip nodes
@@ -179,14 +166,14 @@ let calcBBox nodes =
        maxX = Array.max x
        minY = Array.min y
        maxY = Array.max y
-        center = 0f, 0f
+        center = 0, 0
    }
    |> fun x -> { x with center = calcBBoxCenter x }

 let printBBox grid = calcBBox grid.Nodes |> printfn "%A"

 let bboxToLngLat (coordsys: CoordinateSystem) b =
-    let toLatLon = makeTransform coordsys Transformations.WGS84
+    let toLatLon = makeTransform coordsys CoordSys.WGS84
    let x0, y0 = toLatLon.project ((b.minX, b.minY))
    let x1, y1 = toLatLon.project ((b.maxX, b.maxY))
    // Log.Error (sprintf "(%f, %f) (%f, %f)" x0 y0 x1 y1)
@@ -201,45 +188,68 @@ let bboxToLngLat (coordsys: CoordinateSystem) b =
 let projectBBox proj b =
    let x0, y0 = proj (b.minX, b.minY)
    let x1, y1 = proj (b.maxX, b.maxY)
-    { minX = x0; maxX = x1; minY = y0; maxY = y1; center = proj b.center }
-
-let projectGrid proj (grid: Grid) : Grid =
-    { grid with
-        Nodes = grid.Nodes |> Array.Parallel.map proj
-        BBox = projectBBox proj grid.BBox
+    {
+        minX = x0
+        maxX = x1
+        minY = y0
+        maxY = y1
+        center = proj b.center
    }

+let projectGrid proj (grid: Grid) : Grid = {
+    grid with
+        Nodes = grid.Nodes |> Array.Parallel.map proj
+        BBox = projectBBox proj grid.BBox
+}
+
+let rescaleGrid (factor: float) (grid: Grid) : Grid =
+    let nodes' = grid.Nodes |> Array.Parallel.map (fun (x, y) -> factor * x, factor * y)
+    let bbox' = calcBBox nodes'
+    { grid with Nodes = nodes'; BBox = bbox' }
+
+let translateGrid (x0, y0) grid =
+    let nodes' = grid.Nodes |> Array.Parallel.map (fun (x, y) -> x + x0, y + y0)
+    let bbox' = calcBBox nodes'
+    { grid with Nodes = nodes'; BBox = bbox' }
+
 let toWebMercator (coordsys: CoordinateSystem) (grid: Grid) =
-    let toWebMercator = makeTransform coordsys Transformations.EPSG3857
-    let s = System.Diagnostics.Stopwatch.StartNew()
+    let toWebMercator = makeTransform coordsys CoordSys.EPSG3857
+    let s = System.Diagnostics.Stopwatch.StartNew ()
    let g = grid |> projectGrid toWebMercator.project
-    s.Stop()
+    s.Stop ()
+    Log.Debug $"Reprojected grid: {s.ElapsedMilliseconds} ms"
+    g
+
+let toLonLat (coordsys: CoordinateSystem) (grid: Grid) =
+    let toLonLat = makeTransform coordsys CoordSys.WGS84
+    let s = System.Diagnostics.Stopwatch.StartNew ()
+    let g = grid |> projectGrid toLonLat.project
+    s.Stop ()
    Log.Debug $"Reprojected grid: {s.ElapsedMilliseconds} ms"
    g

 let private chomp (l: string) = l.Split ' ' |> Array.filter ((<>) "")

 // TODO: pattern match and compare sizes for more flex
-let private readGrdHeader (h: string []) =
+let private readGrdHeader (h: string[]) =
    let parse x = (chomp x)[3] |> int
    try
        let nodes = parse h[0]
        let nele = parse h[1]
-        Some(nodes, nele)
-    with
-    | _ -> None
+        Some (nodes, nele)
+    with _ ->
+        None

 let private readObcHeader (h: string) =
    try
        (chomp h)[4] |> int |> Some
-    with
-    | _ -> None
+    with _ ->
+        None

-let private reader (parser: string [] -> 'a) (f: string array) =
+let private reader (parser: string[] -> 'a) (f: string array) =
    try
        f |> Array.map (chomp >> parser) |> Some
-    with
-    | e ->
+    with e ->
        Log.Error e.Message
        None

@@ -254,17 +264,17 @@ let private readElem (f: string array) =
 let private readNodes (f: string array) =
    let p =
        function
-        | [| _; x; y; _ |] -> single x, single y
-        | [| x; y |] -> single x, single y
+        | [| _; x; y; _ |] -> float x, float y
+        | [| x; y |] -> float x, float y
        | x -> failwith $"readNodes failed: %A{x}"
    reader p f

 let private readObc (f: string array) =
    let p =
        function
-        | [| _; x; _ |] -> int x
-        | [| x; _ |] -> int x
-        | [| x |] -> int x
+        | [| _; x; _ |] -> int x - 1
+        | [| x; _ |] -> int x - 1
+        | [| x |] -> int x - 1
        | _ -> failwith "readObc failed"
    reader p f

@@ -276,15 +286,14 @@ let readGrdFile (filename: string) =
        let els, nds = Array.splitAt nele rest
        let elem = readElem els
        let nodes = readNodes nds
-        let toGrid e n =
-            { Grid.empty with Elem = e; Nodes = n; BBox = calcBBox n }
-        toGrid <!> elem <*> nodes)
+        let toGrid e n = { Elem = e; Nodes = n; BBox = calcBBox n }
+        toGrid <!> elem <*> nodes
+    )

 let readObcFile (filename: string) =
    let f = System.IO.File.ReadAllLines filename
    let hdr, rest = Array.splitAt 1 f
-    readObcHeader hdr[0]
-    |> Option.bind (fun _ -> readObc rest)
+    readObcHeader hdr[0] |> Option.bind (fun _ -> readObc rest)

 module Boundary =
    let normalizeElement (a, b, c) =
@@ -298,11 +307,8 @@ module Boundary =
            match Map.tryFind edge a with
            | Some v -> Map.add edge (n :: v) a
            | None -> Map.add edge [ n ] a
-        let appendEdges a (n, x: Edge []) =
-            a
-            |> appendEdge (n, x[0])
-            |> appendEdge (n, x[1])
-            |> appendEdge (n, x[2])
+        let appendEdges a (n, x: Edge[]) =
+            a |> appendEdge (n, x[0]) |> appendEdge (n, x[1]) |> appendEdge (n, x[2])
        grid.Elem
        |> Array.mapi normElIdx
        |> Array.fold appendEdges Map.empty
@@ -320,8 +326,7 @@ module Boundary =
        |> Map.mapValues (fun x -> List.head x, List.last x)

    let makeBoundaryByElementMap (edgeMap: Map<Edge, ElemIdx>) =
-        edgeMap
-        |> Map.fold (fun a k v -> Map.add v k a) Map.empty
+        edgeMap |> Map.fold (fun a k v -> Map.add v k a) Map.empty

    let getBoundaryNodesArray (edgeMap: Map<Edge, ElemIdx>) =
        edgeMap
@@ -331,7 +336,7 @@ module Boundary =
        |> Array.distinct

 module Util =
-    let private floatify = bimap float float
+    // let private floatify = bimap float float

    type Element =
        static member calcCentroid((x0, y0), (x1, y1), (x2, y2)) =
@@ -341,32 +346,65 @@ module Util =
            (x0 + x1 + x2) / 3f, (y0 + y1 + y2) / 3f

        static member calcArea((x0, y0), (x1, y1), (x2, y2)) =
-            x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
-            |> (*) 0.5
+            x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1) |> (*) 0.5

-        static member calcArea((x0, y0), (x1, y1), (x2, y2)) =
-            x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
-            |> (*) 0.5f
+        // static member calcArea((x0, y0), (x1, y1), (x2, y2)) =
+        //     x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
+        //     |> (*) 0.5f

-        static member propToNodal (nIdx: NeighborIndex) (s: single []) =
-            [| 0 .. nIdx.NodesAroundNode.Count - 1 |]
+        static member calcCircumscribedCircle((Ax, Ay as A), (Bx, By as B), (Cx, Cy as C)) =
+            let square (i, j) = i * i + j * j
+            let det (i0, i1, i2) (j0, j1, j2) (k0, k1, k2) =
+                i0 * (j1 * k2 - k1 * j2) - j0 * (i1 * k2 - k1 * i2) + k0 * (i1 * j2 - j1 * i2)
+            let T2 = square A, square B, square C
+            let Tx = Ax, Bx, Cx
+            let Ty = Ay, By, Cy
+            let E = 1.0, 1.0, 1.0
+            let Sx = det T2 Ty E |> (*) 0.5
+            let Sy = det Tx T2 E |> (*) 0.5
+            let a = det Tx Ty E
+            let b = det Tx Ty T2
+            if abs a < 1.0e-12 then
+                failwith "co-linear vertices"
+            let S2 = (Sx, Sy) |> square
+            let r = (b / a + S2 / (a * a)) |> sqrt
+            Sx / a, Sy / a, r
+
+        static member propToNodal (nIdx: NeighborIndex) (s: float[]) =
+            [| 0 .. nIdx.NodesAroundNode.Count - 1 |] // total number of nodes
            |> Array.Parallel.map (fun i ->
                let ns = nIdx.ElemsAroundNode[i]
-                ns
-                |> Array.fold (fun a n -> s[n] + a) 0f
-                |> (*) (1f / single ns.Length))
+                ns |> Array.fold (fun a n -> s[n] + a) 0. |> (*) (1. / float ns.Length)
+            )

-        static member speedToNodal (nIdx: NeighborIndex) (s: (single * single) []) =
-            [| 0 .. nIdx.NodesAroundNode.Count - 1 |]
+        static member speedToNodal (nIdx: NeighborIndex) (s: (float * float)[]) =
+            [| 0 .. nIdx.NodesAroundNode.Count - 1 |] // total number of nodes
            |> Array.Parallel.map (fun i ->
                let ns = nIdx.ElemsAroundNode[i]
                ns
                |> Array.fold
                    (fun a n ->
                        let u, v = s[n]
-                        sqrt (u * u + v * v) + a)
-                    0f
-                |> (*) (1f / single ns.Length))
+                        sqrt (u * u + v * v) + a
+                    )
+                    0.
+                |> (*) (1. / float ns.Length)
+            )
+
+        static member velocityToNodal (nIdx: NeighborIndex) (s: (float * float)[]) =
+            [| 0 .. nIdx.NodesAroundNode.Count - 1 |] // total number of nodes
+            |> Array.Parallel.map (fun i ->
+                let ns = nIdx.ElemsAroundNode[i]
+                let n = float ns.Length
+                ns
+                |> Array.fold
+                    (fun (au, av) n ->
+                        let u, v = s[n]
+                        u + au, v + av
+                    )
+                    (0., 0.)
+                |> fun (u, v) -> u / n, v / n
+            )

    type Node =
        static member calcNodeControlArea (idx: NeighborIndex) (grid: IGrid) =
@@ -377,17 +415,18 @@ module Util =
                getElemsSurroundingNode idx n
                |> Array.map (fun e ->
                    let p0, p1, p2 = grid.getCellVertices e
-                    let x0, y0 as p0' = floatify p0
-                    let x1, y1 = floatify p1
-                    let x2, y2 = floatify p2
+                    let x0, y0 as p0' = p0
+                    let x1, y1 = p1
+                    let x2, y2 = p2
                    let p1' = x0 + half (x1 - x0), y0 + half (y1 - y0)
                    let p2' = x0 + half (x2 - x0), y0 + half (y2 - y0)
-                    let centroid = Element.calcCentroid (p0, p1, p2) |> floatify
+                    let centroid = Element.calcCentroid (p0, p1, p2)
                    let a1 = Element.calcArea (p0', p1', p2')
                    let a2 = Element.calcArea (p1', centroid, p2')
-                    a1 + a2)
+                    a1 + a2
+                )
                |> Array.sum
-                |> single)
+            )

        static member calcNodeArea (idx: NeighborIndex) (grid: IGrid) =
            let nodes = grid.getVertices ()
@@ -395,7 +434,8 @@ module Util =
            |> Array.Parallel.map (fun n ->
                getElemsSurroundingNode idx n
                |> Array.map (grid.getCellVertices >> Element.calcArea)
-                |> Array.sum)
+                |> Array.sum
+            )

    let calcCentroids (grid: IGrid) =
        let n = grid.getVertices ()
@@ -404,43 +444,28 @@ module Util =
            let p0 = n[a]
            let p1 = n[b]
            let p2 = n[c]
-            Element.calcCentroid (p0, p1, p2))
+            Element.calcCentroid (p0, p1, p2)
+        )

    let inline isInsideTriangle (x, y, z) p =
        let sign (p1x, p1y) (p2x, p2y) (p3x, p3y) =
-            (p1x - p3x) * (p2y - p3y)
-            - (p2x - p3x) * (p1y - p3y)
+            (p1x - p3x) * (p2y - p3y) - (p2x - p3x) * (p1y - p3y)

        let d1 = sign p x y
        let d2 = sign p y z
        let d3 = sign p z x

-        let neg = (d1 <= 0.f) || (d2 <= 0.f) || (d3 <= 0.f)
-        let pos = (d1 >= 0.f) || (d2 >= 0.f) || (d3 >= 0.f)
+        let neg = (d1 < 0.) || (d2 < 0.) || (d3 < 0.)
+        let pos = (d1 > 0.) || (d2 > 0.) || (d3 > 0.)

-        neg && pos |> not
-
-    let inline isInsideTriangle_dbl (x, y, z) p =
-        let sign (p1x, p1y) (p2x, p2y) (p3x, p3y) =
-            (p1x - p3x) * (p2y - p3y)
-            - (p2x - p3x) * (p1y - p3y)
-
-        let d1 = sign p x y
-        let d2 = sign p y z
-        let d3 = sign p z x
-
-        let neg = (d1 <= 0.0) || (d2 <= 0.0) || (d3 <= 0.0)
-        let pos = (d1 >= 0.0) || (d2 >= 0.0) || (d3 >= 0.0)
-
-        neg && pos |> not
+        (neg && pos) |> not

    // make a kd-tree for looking up nearest node
    let buildNearestNodeTree (grid: IGrid) =
        grid.getVertices ()
        |> Array.mapi (fun i v ->
            let x, y = v
-            { Pos = { X = x; Y = y; }
-              Data = i }
+            { Pos = x, y; Data = i }
        )
        // |> create2DTree treeLeafSize
        |> createTree
@@ -449,12 +474,10 @@ module Util =
    let buildNearestElementTree (grid: IGrid) =
        grid.getCells ()
        |> Array.mapi (fun i _ ->
-            let pos =
-                i
-                |> grid.getCellVertices
-                |> Element.calcCentroid
-                |> fun (x, y) -> { X = x; Y = y }
-            { Pos = pos; Data = i })
+            let pos = i |> grid.getCellVertices |> Element.calcCentroid
+            // |> fun (x, y) -> { X = x; Y = y }
+            { Pos = pos; Data = i }
+        )
        // |> create2DTree treeLeafSize
        |> createTree

@@ -469,28 +492,39 @@ module Util =
    //         then Some leaf.Data
    //         else None)

-    let tryFindElement
-            (grid: IGrid)
-            (tree: KdTree<single, int>)
-            ((p0, p1): single * single) =
-        tree.GetNearestNeighbours([|p0; p1|], 1)
+    let tryFindElement (grid: IGrid) (tree: KdTree<float, int>) ((p0, p1): float * float) =
+        tree.GetNearestNeighbours ([| p0; p1 |], 1)
        |> Array.tryHead
        |> Option.bind (fun leaf ->
            let vx = grid.getCellVertices leaf.Value
-            if isInsideTriangle vx (p0, p1)
-            then Some leaf.Value
-            else None)
+            if isInsideTriangle vx (p0, p1) then
+                Some leaf.Value
+            else
+                None
+        )

 // type private IdxTree = Tree<Leaf<single, ElemIdx> array, Node<single>>
 // type private NodeIdxTree = Tree<Leaf<single, NodeIdx> array, Node<single>>
-type private IdxTree = KdTree<single, int>
-type private NodeIdxTree = KdTree<single, int>
+type private IdxTree = KdTree<float, int>
+type private NodeIdxTree = KdTree<float, int>
+
+[<MessagePackObject>]
+type BinGrid = {
+    [<Key(0)>]
+    hash: byte[]
+    [<Key(1)>]
+    vertices: (float * float)[]
+    [<Key(2)>]
+    cells: (int * int * int)[]
+} with
+    member this.toGrid() : Grid = { Nodes = this.vertices; Elem = this.cells; BBox = calcBBox this.vertices }

 type ExtendedGrid(grid: IGrid) =
-    let mutable nodeTree : NodeIdxTree option = None
-    let mutable elementTree : IdxTree option = None
-    let mutable neighborIndex : NeighborIndex option = None
-    let mutable centroids : Vertex [] option = None
+    let mutable nodeTree: NodeIdxTree option = None
+    let mutable elementTree: IdxTree option = None
+    let mutable neighborIndex: NeighborIndex option = None
+    let mutable centroids: Vertex[] option = None
+    let mutable gridHash: byte[] = [||]

    let getNeighborIdx () =
        match neighborIndex with
@@ -500,97 +534,117 @@ type ExtendedGrid(grid: IGrid) =
            neighborIndex.Value

    interface IGrid with
-        member x.getVertex n = grid.getVertex n
-        member x.getCell n = grid.getCell n
-        member x.getCellVertices n = grid.getCellVertices n
-        member x.getVertices() = grid.getVertices ()
-        member x.getCells() = grid.getCells ()
-        member x.getBoundingBox() = grid.getBoundingBox ()
-        member x.Projection = grid.Projection
+        member this.getVertex n = grid.getVertex n
+        member this.getCell n = grid.getCell n
+        member this.getCellVertices n = grid.getCellVertices n
+        member this.getVertices() = grid.getVertices ()
+        member this.getCells() = grid.getCells ()
+        member this.getBoundingBox() = grid.getBoundingBox ()

-    member x.NeighborIndex = neighborIndex
-    member x.NodeTree = nodeTree
-    member x.ElementTree = elementTree
+    member this.NeighborIndex = neighborIndex
+    member this.NodeTree = nodeTree
+    member this.ElementTree = elementTree

-    member x.Grid = grid
+    member this.Grid = grid

-    member x.initNeighborIndex(?cache: string) =
+    member this.ToGrid() = {
+        Nodes = this.Grid.getVertices ()
+        Elem = this.Grid.getCells ()
+        BBox = this.Grid.getBoundingBox ()
+    }
+
+
+    member this.initNeighborIndex(?cache: string) =
        match neighborIndex with
        | Some _ -> ()
        | None -> neighborIndex <- makeNeighborIndex grid |> Some

-    member x.initNodeTree() =
+    member this.initNodeTree() =
        match nodeTree with
        | Some _ -> ()
        | None -> nodeTree <- Util.buildNearestNodeTree grid |> Some

-    member x.initElementTree() =
+    member this.initElementTree() =
        match elementTree with
        | Some _ -> ()
        | None -> elementTree <- Util.buildNearestElementTree grid |> Some

-    // member x.nearestNode(p0: single, p1: single) =
-    //     let nearest (tree: Tree<_, _>) =
-    //         nearestNeighbor tree { X = p0; Y = p1 }
-    //         |> Option.map (fun l -> l.Data)
-    //     match nodeTree with
-    //     | Some tree -> nearest tree
-    //     | None ->
-    //         x.initNodeTree ()
-    //         nearest nodeTree.Value
-
-    member x.nearestNode(p0: single, p1: single) =
+    member this.nearestNode(p0: float, p1: float) =
        let nearest (tree: KdTree<_, _>) =
-            tree.GetNearestNeighbours([|p0; p1|], 1)
+            tree.GetNearestNeighbours ([| p0; p1 |], 1)
            |> Array.tryHead
            |> Option.map (fun l -> l.Value)
        match nodeTree with
        | Some tree -> nearest tree
        | None ->
-            x.initNodeTree ()
+            this.initNodeTree ()
            nearest nodeTree.Value

-    // member private x.tryFindElementTwice (grid: IGrid) (tree: IdxTree) ((p0, p1): single * single as p) =
-    //     match Util.tryFindElement grid tree p with
-    //     | Some elemIdx -> Some elemIdx
-    //     | None ->
-    //         nearestNeighbor tree { X = p0; Y = p1 }
-    //         |> Option.bind (fun leaf ->
-    //             x.getElemsSurroundingElem leaf.Data
-    //             |> Array.tryFind (fun eIdx ->
-    //                 let vx = x.Grid.getCellVertices eIdx
-    //                 Util.isInsideTriangle vx p))
+    member this.tryGetNode(p: float * float) =
+        this.nearestNode p
+        |> Option.bind (fun n ->
+            this.getElemsSurroundingNode n
+            |> Array.fold
+                (fun (a: int option) e ->
+                    if a.IsNone then
+                        let vx = this.Grid.getCellVertices e
+                        if Util.isInsideTriangle vx p then Some n else None
+                    else
+                        a
+                )
+                None
+        )

-    member private x.tryFindElementTwice (grid: IGrid) (tree: IdxTree) ((p0, p1): single * single as p) =
-        match Util.tryFindElement grid tree p with
-        | Some elemIdx -> Some elemIdx
-        | None ->
-            tree.GetNearestNeighbours([|p0; p1|], 1)
+    member private this.tryFindElementTwice (grid: IGrid) (tree: IdxTree) ((p0, p1): float * float as p) =
+        Util.tryFindElement grid tree p
+        |> Option.orElse (
+            tree.GetNearestNeighbours ([| p0; p1 |], 1)
            |> Array.tryHead
            |> Option.bind (fun leaf ->
-                x.getElemsSurroundingElem leaf.Value
+                this.getElemsSurroundingElem leaf.Value
                |> Array.tryFind (fun eIdx ->
-                    let vx = x.Grid.getCellVertices eIdx
-                    Util.isInsideTriangle vx p))
+                    let vx = this.Grid.getCellVertices eIdx
+                    Util.isInsideTriangle vx p
+                )
+            )
+        )

-    member x.tryGetElement p =
-        match elementTree with
-        | Some tree -> x.tryFindElementTwice grid tree p
-        | None ->
-            x.initElementTree ()
-            Log.Debug $"You shouldn't be here!"
-            // x.initNodeTree ()
-            x.tryFindElementTwice grid elementTree.Value p
+    member this.tryGetElement p =
+        elementTree
+        |> Option.bind (fun tree -> this.tryFindElementTwice grid tree p)
+        |> Option.orElseWith (fun () ->
+            this.initElementTree ()
+            this.tryFindElementTwice grid elementTree.Value p
+        )

-    member x.tryGetElementSloppy(p0, p1 as p) =
+    member this.tryGetElementSloppy(p0, p1 as p) =
        match elementTree with
        | Some tree -> Util.tryFindElement grid tree p
        | None ->
-            x.initElementTree ()
-            // x.initNodeTree ()
+            this.initElementTree ()
+            // this.initNodeTree ()
            Util.tryFindElement grid elementTree.Value p

-    member x.getCentroids() =
+    member this.initHash(hash: byte[]) =
+        if gridHash.Length = 0 && hash.Length > 0 then
+            gridHash <- hash
+
+    member this.getHash() =
+        if gridHash.Length = 0 then
+            let bg: BinGrid = {
+                hash = [||]
+                vertices = (this :> IGrid).getVertices ()
+                cells = (this :> IGrid).getCells ()
+            }
+            let bytes = MessagePackSerializer.Serialize (bg)
+            let sha1 = System.Security.Cryptography.SHA1.Create ()
+            gridHash <- sha1.ComputeHash bytes
+        gridHash
+
+    member this.getHashString() =
+        this.getHash () |> Convert.ToHexStringLower
+
+    member this.getCentroids() =
        match centroids with
        | Some cx -> cx
        | None ->
@@ -598,54 +652,55 @@ type ExtendedGrid(grid: IGrid) =
            centroids <- Some cx
            cx

-    member x.getElemsSurroundingNode n =
-        getNeighborIdx ()
-        |> fun idx -> idx.ElemsAroundNode[n]
+    member this.calcCircumCircle e =
+        let triangle = grid.getCellVertices e
+        Util.Element.calcCircumscribedCircle triangle

-    member x.getNodesSurroundingNode n =
-        getNeighborIdx ()
-        |> fun idx -> idx.NodesAroundNode[n]
+    member this.getElemsSurroundingNode n =
+        getNeighborIdx () |> fun idx -> idx.ElemsAroundNode[n]

-    member x.getNodesSurroundingElem e =
+    member this.getNodesSurroundingNode n =
+        getNeighborIdx () |> fun idx -> idx.NodesAroundNode[n]
+
+    member this.getNodesSurroundingElem e =
        let idx = getNeighborIdx ()
        let elem = grid.getCell e
        getSurrounding idx.NodesAroundNode elem

-    member x.getElemsSurroundingElem e =
+    member this.getElemsSurroundingElem e =
        let idx = getNeighborIdx ()
        let elem = grid.getCell e
        getSurrounding idx.ElemsAroundNode elem

-    member x.saveNeighborIndex(fname: string) =
-        let binarySerializer = FsPickler.CreateBinarySerializer()
+    member this.saveNeighborIndex(fname: string) =
+        let binarySerializer = FsPickler.CreateBinarySerializer ()
        let nix = getNeighborIdx ()
        let pickle = binarySerializer.Pickle nix
-        IO.File.WriteAllBytes(fname, pickle)
+        IO.File.WriteAllBytes (fname, pickle)

-    member x.loadNeighborIndex(fname: string) =
-        let binarySerializer = FsPickler.CreateBinarySerializer()
+    member this.loadNeighborIndex(fname: string) =
+        let binarySerializer = FsPickler.CreateBinarySerializer ()
        if IO.File.Exists fname then
            let pickle = IO.File.ReadAllBytes fname
-            neighborIndex <-
-                binarySerializer.UnPickle<NeighborIndex> pickle
-                |> Some
+            neighborIndex <- binarySerializer.UnPickle<NeighborIndex> pickle |> Some
            true
        else
            false

-    member x.saveNodeTree(fname: string) =
-        let binarySerializer = FsPickler.CreateBinarySerializer()
+    member this.saveNodeTree(fname: string) =
+        let binarySerializer = FsPickler.CreateBinarySerializer ()
        let tree =
            nodeTree
            |> Option.defaultWith (fun () ->
-                x.initNodeTree ()
-                nodeTree.Value)
+                this.initNodeTree ()
+                nodeTree.Value
+            )

        let pickle = binarySerializer.Pickle tree
-        IO.File.WriteAllBytes(fname, pickle)
+        IO.File.WriteAllBytes (fname, pickle)

-    member x.loadNodeTree(fname: string) =
-        let binarySerializer = FsPickler.CreateBinarySerializer()
+    member this.loadNodeTree(fname: string) =
+        let binarySerializer = FsPickler.CreateBinarySerializer ()
        if IO.File.Exists fname then
            let pickle = IO.File.ReadAllBytes fname
            nodeTree <- binarySerializer.UnPickle<IdxTree> pickle |> Some
@@ -653,22 +708,23 @@ type ExtendedGrid(grid: IGrid) =
        else
            false

-    member x.saveElementTree(fname: string) =
-        let binarySerializer = FsPickler.CreateBinarySerializer()
+    member this.saveElementTree(fname: string) =
+        let binarySerializer = FsPickler.CreateBinarySerializer ()
        let tree =
            elementTree
            |> Option.defaultWith (fun () ->
-                x.initElementTree ()
-                elementTree.Value)
+                this.initElementTree ()
+                elementTree.Value
+            )

        let pickle = binarySerializer.Pickle tree
        IO.File.WriteAllBytes (fname, pickle)

-    member x.loadElementTree(fname: string) =
-        let binarySerializer = FsPickler.CreateBinarySerializer()
+    member this.loadElementTree(fname: string) =
+        let binarySerializer = FsPickler.CreateBinarySerializer ()
        if IO.File.Exists fname then
            let pickle = IO.File.ReadAllBytes fname
            elementTree <- binarySerializer.UnPickle<IdxTree> pickle |> Some
            true
        else
-            false
+            false
--- a/src/Interpol.fs
+++ b/src/Interpol.fs
@@ -8,8 +8,8 @@ type InterpolCoefs = ((int * int) * (float * float)) array array

 type DepthInterpolCoefs = { iRho: InterpolCoefs; iU: InterpolCoefs; iV: InterpolCoefs }

-let private findNearestZ z0 (h: float []) =
-    let rec findNerest' z0 (h: float []) d n =
+let private findNearestZ z0 (h: float[]) =
+    let rec findNerest' z0 (h: float[]) d n =
        if n < h.Length - 1 then
            let d' = abs (h[n] - z0)
            if d' < d then findNerest' z0 h d' (n + 1) else n
@@ -21,9 +21,10 @@ let private calcInterpolationWeightedIdx (rn, fn) =
    let nearestIdx = Array.map (fun fzi -> findNearestZ fzi rn) fn
    nearestIdx
    |> Array.mapi (fun i j ->
+        let last = Array.last rn
        if fn[i] <= rn[0] then
            (0, 0), (1.0, 0.0)
-        elif fn[i] >= rn[^0] then
+        elif fn[i] >= last then
            (rn.Length - 1, 0), (1.0, 0.0)
        elif abs (fn[i] - rn[j]) < 9.999999975e-07 then
            (j, 0), (1.0, 0.0)
@@ -40,42 +41,36 @@ let private calcInterpolationWeightedIdx (rn, fn) =
        else
            failwith "not reachable")

-let private calcInterpMatrices (rz: float [] []) (fz: float [] []) =
-    Array.zip rz fz
-    |> Array.Parallel.map calcInterpolationWeightedIdx
+let private calcInterpMatrices (rz: float[][]) (fz: float[][]) =
+    Array.zip rz fz |> Array.Parallel.map calcInterpolationWeightedIdx

-let mkDepthInterpolCoefs (siglay: single [,]) (h: float []) (roms: float [] []) =
-    let s =
-        siglay
-        |> Array2D.map float
-        |> Matrix.Build.DenseOfArray
+let mkDepthInterpolCoefs (siglay: single[,]) (h: float[]) (roms: float[][]) =
+    let s = siglay |> Array2D.map float |> Matrix.Build.DenseOfArray
    let rescale (sigma: float Matrix) h =
        let sh = h |> Array.map float
-        sigma.MapIndexed(fun _ j x -> -x * sh[j])
-        |> fun x -> x.ToColumnArrays()
+        sigma.MapIndexed (fun _ j x -> -x * sh[j]) |> fun x -> x.ToColumnArrays ()

    let conv m =
        m
        |> Array.map (Array.map ((*) -1.0))
        |> matrix
-        |> fun x -> x.ToColumnArrays()
+        |> fun x -> x.ToColumnArrays ()

    let romz = conv roms
    let z = rescale s h
    calcInterpMatrices romz z

-let zInterpolProp (iz: InterpolCoefs) (adjRomsProp: float [] []) =
+let zInterpolProp (iz: InterpolCoefs) (adjRomsProp: float[][]) =
    let pz = adjRomsProp |> matrix
    iz
    |> Array.Parallel.mapi (fun n x ->
-        let p = pz[ *, n ].ToArray() |> Array.rev
-        x
-        |> Array.map (fun ((i1, i2), (w1, w2)) -> p[i1] * w1 + p[i2] * w2))
+        let p = pz[*, n].ToArray () |> Array.rev
+        x |> Array.map (fun ((i1, i2), (w1, w2)) -> p[i1] * w1 + p[i2] * w2))
    |> matrix
-    |> fun x -> x.Transpose().ToArray()
+    |> fun x -> x.Transpose().ToArray ()

 // Trivially already satisfied by ajoint
-let hInterpolNearestProp (adjRomsProp: float [] []) = adjRomsProp
+let hInterpolNearestProp (adjRomsProp: float[][]) = adjRomsProp

 type BiWght = float * float * float * float
 type BiProp = float * float * float * float
@@ -88,51 +83,47 @@ let genBilinearInterpolationWgts
    : BiWght option =
    if mask0 && mask1 && mask2 && mask3 then
        let AB = (x1 - x0) ** 2.0 + (y1 - y0) ** 2.0 |> sqrt
-        let BC = (x2 - x1)**2.0 + (y2 - y1)**2.0 |> sqrt
+        let BC = (x2 - x1) ** 2.0 + (y2 - y1) ** 2.0 |> sqrt
        let CD = (x3 - x2) ** 2.0 + (y3 - y2) ** 2.0 |> sqrt
        let DA = (x0 - x3) ** 2.0 + (y0 - y3) ** 2.0 |> sqrt
        let h0 =
-            let a0 =
-                (x * (y0 - y3) + x0 * (y3 - y) + x3 * (y - y0))
-                * 0.5
+            let a0 = (x * (y0 - y3) + x0 * (y3 - y) + x3 * (y - y0)) * 0.5
            2.0 * a0 / DA
        let h1 =
-            let a1 =
-                (x * (y1 - y0) + x0 * (y - y1) + x1 * (y0 - y))
-                * 0.5
+            let a1 = (x * (y1 - y0) + x0 * (y - y1) + x1 * (y0 - y)) * 0.5
            2.0 * a1 / AB
        let w00 = (h0 / AB - 1.0) * (h1 / DA - 1.0)
        let w10 = h0 / AB * (1.0 - h1 / DA)
        let w11 = h1 / DA * h0 / CD
        let w01 = h1 / DA * (1.0 - h0 / CD)
-        Some(w00, w10, w11, w01)
+        Some (w00, w10, w11, w01)
    else
        None

 let makeBiWeights
-    (pos: (float * float) [])
-    (boxes: ((float * float) * (float * float) * (float * float) * (float * float)) [])
-    (mask: (bool * bool * bool * bool) [])
+    (pos: (float * float)[])
+    (boxes: ((float * float) * (float * float) * (float * float) * (float * float))[])
+    (mask: (bool * bool * bool * bool)[])
    =
    Array.zip3 pos boxes mask
    |> Array.map (fun (p, b, m) -> genBilinearInterpolationWgts p b m)

-let interpolateCells (wghts: BiWght option []) (prop: BiProp []) =
+let interpolateCells (wghts: BiWght option[]) (prop: BiProp[]) =
    Array.zip wghts prop
    |> Array.Parallel.map (fun (w, (p00, p01, p11, p10)) ->
        match w with
        | Some w ->
            let w00, w01, w11, w10 = w
-            Some(w00 * p00 + w01 * p01 + w11 * p11 + w10 * p10)
+            Some (w00 * p00 + w01 * p01 + w11 * p11 + w10 * p10)
        | None -> None)
 //    |> Array.Parallel.map (fun ((w00, w01, w11, w10), (p00, p01, p11, p10)) ->
 //        w00 * p00 + w01 * p01 + w11 * p11 + w10 * p10)

-let private valtest (intval: float option []) =
+let private valtest (intval: float option[]) =
    let somval = intval |> Array.filter (fun v -> Option.isSome v)
    somval.Length = intval.Length

-let private fillOutOfBounds (grid: Grid) (interpVal: float option []) =
+let private fillOutOfBounds (grid: Grid) (interpVal: float option[]) =
    let node = interpVal.Length = grid.Nodes.Length
    let nbridx = makeNeighborIndex (grid :> IGrid)
    let mutable allval = valtest interpVal
@@ -141,24 +132,20 @@ let private fillOutOfBounds (grid: Grid) (interpVal: float option []) =
            if Option.isNone interpVal[n] then
                let nb =
                    if node then
-                        nbridx.NodesAroundNode[n]
-                        |> Array.filter (fun k -> Option.isSome interpVal[k])
+                        nbridx.NodesAroundNode[n] |> Array.filter (fun k -> Option.isSome interpVal[k])
                    else
                        getElemsSurroundingElem nbridx grid n
                        |> Array.filter (fun k -> Option.isSome interpVal[k])
                if nb.Length > 0 then
-                    let nval =
-                        nb
-                        |> Array.map (fun k -> interpVal[k].Value)
-                        |> Array.average
+                    let nval = nb |> Array.map (fun k -> interpVal[k].Value) |> Array.average
                    interpVal[n] <- Some nval
        allval <- valtest interpVal
    interpVal |> Array.map (fun v -> v.Value)

 let interpolateProp
    (coords, _ as grid: BiPos * Mask)
-    (pos: (float * float) [])
-    (prop: float [,])
+    (pos: (float * float)[])
+    (prop: float[,])
    (fvgrid: Grid)
    (oob: bool)
    =
@@ -179,8 +166,8 @@ let interpolateProp
            iprop
            |> Array.mapi (fun id p ->
                match p with
-                |Some _ -> -1
-                |None -> id)
+                | Some _ -> -1
+                | None -> id)
            |> Array.filter (fun i -> i > 0)
            |> Array.fold (fun a i -> a + "; " + string i) ""
        failwith ($"FVCOM grid out of bounds. Oob indexes: {sidx}")
@@ -189,7 +176,7 @@ let interpolateProp
    else
        iprop |> Array.map (fun v -> v.Value)

-let interpCoefs (coords, _ as grid: BiPos * Mask) (pos: (float * float) []) =
+let interpCoefs (coords, _ as grid: BiPos * Mask) (pos: (float * float)[]) =
    let cm = getNearestCellCorner grid pos
    let boxes, mask =
        cm
@@ -198,12 +185,10 @@ let interpCoefs (coords, _ as grid: BiPos * Mask) (pos: (float * float) []) =
        |> Array.unzip
    let wgths = makeBiWeights pos boxes mask
    let corneridx =
-        cm
-        |> Array.map (fun c -> Array.unzip c)
-        |> Array.map (fun (i, _) -> i)
+        cm |> Array.map (fun c -> Array.unzip c) |> Array.map (fun (i, _) -> i)
    corneridx, wgths

-let interpProp (corneridx: (int * int) [] []) (wgths: BiWght option []) (prop: float [,]) (grid: Grid) (oob: bool) =
+let interpProp (corneridx: (int * int)[][]) (wgths: BiWght option[]) (prop: float[,]) (grid: Grid) (oob: bool) =
    let boxProps =
        corneridx
        |> Array.map (fun box -> box |> Array.map (fun (n, m) -> prop[n, m]))
@@ -214,8 +199,8 @@ let interpProp (corneridx: (int * int) [] []) (wgths: BiWght option []) (prop: f
            iprop
            |> Array.mapi (fun i p ->
                match p with
-                |Some _ -> -1
-                |None -> i)
+                | Some _ -> -1
+                | None -> i)
            |> Array.filter (fun i -> i >= 0)
            |> Array.fold (fun a i -> a + "; " + string i) ""
        failwith ($"FVCOM grid out of bounds. Oob indexes: {sidx}")
--- a/src/NorKyst.fs
+++ b/src/NorKyst.fs
@@ -1,18 +1,21 @@
 module Oceanbox.FvcomKit.NorKyst

 open System
+
 open Thredds

 let private getNorkystUrl (threddsUrl: string) (d: DateTime) =
    let url = threddsUrl
    let fmt (d: DateTime) =
        $"{d.Year}%02d{d.Month}%02d{d.Day}T00Z.nc"
+
    $"{url}/fou-hi/new_norkyst800m/his/ocean_his.an.{fmt d}"

 let private getArchiveUrls urlf (t: DateTime) =
-    let t = t.ToUniversalTime()
-    let now = DateTime.Now.ToUniversalTime()
+    let t = t.ToUniversalTime ()
+    let now = DateTime.Now.ToUniversalTime ()
    let dDay = (now.Date - t.Date).Days
+
    if dDay < 0 then // no data available
        []
    else
@@ -21,6 +24,4 @@ let private getArchiveUrls urlf (t: DateTime) =
 let tryGetArchive (threddsUrl: string) (t: DateTime) =
    getArchiveUrls (getNorkystUrl threddsUrl) t
    |> tryOpenThredds
-    |> Option.bind (fun (_, ds) ->
-        tryGetTimeIndex ds t
-        |> Option.map (fun idx -> (ds, idx)))
+    |> Option.bind (fun (_, ds) -> tryGetTimeIndex ds t |> Option.map (fun idx -> ds, idx))
--- a/src/NorShelf.fs
+++ b/src/NorShelf.fs
@@ -10,12 +10,14 @@ let private getNorshelfUrl (threddsUrl: string) (kind: Kind) (mode: Mode) (date:
    let fmt (d: DateTime) =
        $"{d.Year}%02d{d.Month}%02d{d.Day}T00Z.nc"
    let url = threddsUrl + "/sea_norshelf_files"
-    $"{url}/norshelf_{kind}_{mode}_{fmt date}"
+
+    $"{url}/{date.Year}/%02d{date.Month}/norshelf_{kind}_{mode}_{fmt date}"

 let private getArchiveUrls urlf (t: DateTime) =
-    let t = t.ToUniversalTime()
-    let now = DateTime.Now.ToUniversalTime()
+    let t = t.ToUniversalTime ()
+    let now = DateTime.Now.ToUniversalTime ()
    let dDay = (now.Date - t.Date).Days
+
    if dDay < -3 then // no data available
        []
    elif dDay <= 0 then // forecast, count down from latest
@@ -27,11 +29,10 @@ let private getArchiveUrls urlf (t: DateTime) =

 let tryGetArchive threddsUrl avg (t: DateTime) =
    let kind = if avg then Avg else Qck
+
    getArchiveUrls (getNorshelfUrl threddsUrl kind) t
    |> tryOpenThredds
-    |> Option.bind (fun (fname, ds) ->
-        tryGetTimeIndex ds t
-        |> Option.map (fun idx -> (fname, ds, idx)))
+    |> Option.bind (fun (fname, ds) -> tryGetTimeIndex ds t |> Option.map (fun idx -> (fname, ds, idx)))

 let readArchive file reader =
    tryOpenArchive file
@@ -40,8 +41,7 @@ let readArchive file reader =
            let data = reader nc
            Log.Information $"Read NorShelf data from {file}"
            Some data
-        with
-        | e ->
+        with e ->
            Log.Error e.Message
            None)

--- a/src/Oceanbox.FvcomKit.fsproj
+++ b/src/Oceanbox.FvcomKit.fsproj
@@ -2,17 +2,19 @@
 <Project Sdk="Microsoft.NET.Sdk">
  <PropertyGroup>
    <OutputType>Library</OutputType>
-    <TargetFramework>net6.0</TargetFramework>
-    <IsPackable>true</IsPackable>
-    <PackageId>Oceanbox.FvcomKit</PackageId>
+    <TargetFramework>net9.0</TargetFramework>
+    <Company>Oceanbox AS</Company>
    <Authors/>
-    <Company/>
-    <Version>3.12.2</Version>
-    <LangVersion>preview</LangVersion>
+    <PackageId>Oceanbox.FvcomKit</PackageId>
+    <RuntimeIdentifier>linux-x64</RuntimeIdentifier>
+    <Version>6.0.0</Version>
+    <RestorePackagesWithLockFile>false</RestorePackagesWithLockFile>
+    <IsPackable>true</IsPackable>
  </PropertyGroup>
  <ItemGroup>
    <Compile Include="Types.fs"/>
    <Compile Include="Grid.fs"/>
+    <Compile Include="Arome.fs"/>
    <Compile Include="Fvcom.fs"/>
    <Compile Include="Thredds.fs"/>
    <Compile Include="Polygon.fs"/>
@@ -22,24 +24,23 @@
    <Compile Include="NorKyst.fs"/>
    <Compile Include="NorShelf.fs"/>
    <Compile Include="Smoothing.fs"/>
-    <Compile Include="Gradient.fs"/>
-    <Compile Include="Evaluate.fs"/>
  </ItemGroup>
  <ItemGroup>
-    <PackageReference Update="FSharp.Core" Version="6.0.5"/>
-    <PackageReference Include="FSharp.Data" Version="6.0.1-beta001"/>
-    <PackageReference Include="FSharpPlus" Version="1.2.4"/>
+    <PackageReference Include="FSharp.Data" Version="6.4.1"/>
+    <PackageReference Include="FSharpPlus" Version="1.7.0"/>
    <PackageReference Include="FsPickler" Version="5.3.2"/>
    <PackageReference Include="KDTree" Version="1.4.1"/>
    <PackageReference Include="MathNet.Numerics.FSharp" Version="5.0.0"/>
-    <PackageReference Include="sdslite" Version="2.5.0"/>
-    <PackageReference Include="Serilog" Version="2.12.0-dev-01533"/>
-    <PackageReference Include="Serilog.Sinks.Console" Version="4.0.2-dev-00890"/>
-    <PackageReference Include="Serilog.Sinks.Seq" Version="5.1.2-dev-00225"/>
-    <PackageReference Include="Thoth.Json.Net" Version="8.0.0"/>
+    <PackageReference Include="MessagePack" Version="3.1.3"/>
+    <PackageReference Include="ProjNet.FSharp" Version="5.2.0"/>
+    <PackageReference Include="Oceanbox.SDSLite" Version="2.8.0"/>
+    <PackageReference Include="Serilog" Version="4.2.0"/>
+    <PackageReference Include="Serilog.Sinks.Console" Version="6.0.0"/>
+    <PackageReference Include="Serilog.Sinks.Seq" Version="9.0.0"/>
+    <PackageReference Update="FSharp.Core" Version="9.0.303"/>
  </ItemGroup>
  <ItemGroup>
-    <ProjectReference Include="..\submodules\ProjNet.FSharp\src\ProjNet.FSharp.fsproj"/>
-    <ProjectReference Include="..\submodules\FsKdTree\lib\FsKDTree.fsproj"/>
+    <PackageUpdate Include="ProjNet.FSharp" Version="*" Condition=" '$(ContinuousIntegrationBuild)'=='true' "/>
+    <PackageUpdate Include="Oceanbox.SDSLite" Version="*" Condition=" '$(ContinuousIntegrationBuild)'=='true' "/>
  </ItemGroup>
 </Project>
--- a/src/Polygon.fs
+++ b/src/Polygon.fs
@@ -2,7 +2,7 @@

 type Point = float * float

-type Polygon = Point []
+type Polygon = Point[]

 let private maxdist (pol: Polygon) (pnt: Point) =
    pol
@@ -28,10 +28,7 @@ let private onsegment (p: Polygon) =
    let maxy = y[0..1] |> Array.max
    let miny = y[0..1] |> Array.min

-    if x[2] >= minx
-       && x[2] <= maxx
-       && y[2] >= miny
-       && y[2] <= maxy then
+    if x[2] >= minx && x[2] <= maxx && y[2] >= miny && y[2] <= maxy then
        true
    else
        false
@@ -56,8 +53,7 @@ let private intersect (p1: Polygon) (p2: Polygon) =
 let private isodd a = if a % 2 = 0 then false else true

 let inpolygon (pol: Polygon) (point: Point) =
-    let (line: Polygon) =
-        [| fst point, snd point; fst point + (maxdist pol point), snd point |]
+    let (line: Polygon) = [| fst point, snd point; fst point + (maxdist pol point), snd point |]
    let poly =
        if pol[0] <> pol[pol.Length - 1] then
            Array.append pol [| pol[0] |] |> Array.pairwise
--- a/src/ROMS.fs
+++ b/src/ROMS.fs
@@ -11,36 +11,33 @@ open Adjoin
 open Interpol

 type GridData<'a> = { rho: 'a; u: 'a; v: 'a }
-type GridMask = GridData<bool [,]>
+type GridMask = GridData<bool[,]>
 type GridPos = GridData<BiPos>
-type GridVars = GridData<float [,]>
+type GridVars = GridData<float[,]>

-type RomsGrid =
-    {
-        h: GridVars
-        z: GridVars []
-        pos: GridPos
-        wetMask: GridMask
-        angle: float [,]
-    }
+type RomsGrid = {
+    h: GridVars
+    z: GridVars[]
+    pos: GridPos
+    wetMask: GridMask
+    angle: float[,]
+}

-type RomsProps =
-    {
-        salt: float [,] []
-        temp: float [,] []
-        zeta: float [,]
-        u: float [,] []
-        v: float [,] []
-    }
+type RomsProps = {
+    salt: float[,][]
+    temp: float[,][]
+    zeta: float[,]
+    u: float[,][]
+    v: float[,][]
+}

-type AdjoinedProp =
-    {
-        salt: float [] []
-        temp: float [] []
-        zeta: float []
-        u: float [] []
-        v: float [] []
-    }
+type AdjoinedProp = {
+    salt: float[][]
+    temp: float[][]
+    zeta: float[]
+    u: float[][]
+    v: float[][]
+}

 type FvcomGrid = Fvcom.FvcomGrid

@@ -49,58 +46,45 @@ let private layerZ csR (vars: GridVars) =
    let rho = csR |> Array.map (mult vars.rho)
    let u = csR |> Array.map (mult vars.u)
    let v = csR |> Array.map (mult vars.v)
-    Array.zip3 rho u v
-    |> Array.map (fun (a, b, c) -> { rho = a; u = b; v = c })
+    Array.zip3 rho u v |> Array.map (fun (a, b, c) -> { rho = a; u = b; v = c })

 //  convert fractional ROMS landmask to ocean mask
 let readWetMask (nc: DataSet) =
    let isWet = Array2D.map ((=) 1.)
    {
-        rho = nc[ "mask_rho" ].GetData() :?> float [,] |> isWet
-        u = nc[ "mask_u" ].GetData() :?> float [,] |> isWet
-        v = nc[ "mask_v" ].GetData() :?> float [,] |> isWet
+        rho = nc["mask_rho"].GetData () :?> float[,] |> isWet
+        u = nc["mask_u"].GetData () :?> float[,] |> isWet
+        v = nc["mask_v"].GetData () :?> float[,] |> isWet
    }

 let private readProp3 (nc: DataSet) (p: string) t (prop: string) =
    let nEta = nc.Dimensions[$"eta_{p}"].Length
    let nXi = nc.Dimensions[$"xi_{p}"].Length
-    let scaling =
-        nc.Variables[prop].Metadata["scale_factor"] :?> single
-        |> float
-    let offset =
-        nc.Variables[prop].Metadata["add_offset"] :?> single
-        |> float
+    let scaling = nc.Variables[prop].Metadata["scale_factor"] :?> single |> float
+    let offset = nc.Variables[prop].Metadata["add_offset"] :?> single |> float
    let o = [| t; 0; 0 |]
    let shp = [| 1; nEta; nXi |]
-    let convert (x: int16 [,,]) =
-        x[0, *, *]
-        |> Array2D.map (fun x -> float x * scaling + offset)
-    nc[ prop ].GetData(o, shp) :?> int16 [,,]
-    |> convert
+    let convert (x: int16[,,]) =
+        x[0, *, *] |> Array2D.map (fun x -> float x * scaling + offset)
+    nc[prop].GetData (o, shp) :?> int16[,,] |> convert

 let private readProp4 (nc: DataSet) (p: string) t (prop: string) =
    let nEta = nc.Dimensions[$"eta_{p}"].Length
    let nXi = nc.Dimensions[$"xi_{p}"].Length
    let nS = nc.Dimensions["s_rho"].Length
-    let scaling =
-        nc.Variables[prop].Metadata["scale_factor"] :?> single
-        |> float
-    let offset =
-        nc.Variables[prop].Metadata["add_offset"] :?> single
-        |> float
+    let scaling = nc.Variables[prop].Metadata["scale_factor"] :?> single |> float
+    let offset = nc.Variables[prop].Metadata["add_offset"] :?> single |> float
    let o = [| t; 0; 0; 0 |]
    let shp = [| 1; nS; nEta; nXi |]
-    let convert (x: int16 [,,,]) =
+    let convert (x: int16[,,,]) =
        x[0, *, *, *]
        |> Array3D.map (fun x -> float x * scaling + offset)
-        |> fun m ->
-            [|
-                for z = 0 to nS - 1 do
-                    m[z, *, *]
-            |]
+        |> fun m -> [|
+            for z = 0 to nS - 1 do
+                m[z, *, *]
+        |]

-    nc[ prop ].GetData(o, shp) :?> int16 [,,,]
-    |> convert
+    nc[prop].GetData (o, shp) :?> int16[,,,] |> convert

 let readProps (nc: DataSet) t : RomsProps =
    Log.Information "Reading ROMS props..."
@@ -119,32 +103,29 @@ let readVerticalGrid path =
        match Thredds.tryOpenArchive path with
        | Some ds -> ds
        | None -> failwith "open vertical grid failed"
-    nc[ "Cs_r" ].GetData() :?> float [] |> Array.rev
+    nc["s_rho"].GetData () :?> float[] |> Array.rev

 let readGrid (nc: DataSet) =
    Log.Information "Reading ROMS grid..."

-    let hRho =
-        nc[ "h" ].GetData() :?> float [,]
-        |> Matrix.Build.DenseOfArray
+    let hRho = nc["h"].GetData () :?> float[,] |> Matrix.Build.DenseOfArray

-    let csR = nc[ "Cs_r" ].GetData() :?> float [] |> Array.rev
+    let srho = nc["s_rho"].GetData () :?> float[] |> Array.rev

    let e1, e2 = hRho.ColumnCount - 2, hRho.RowCount - 2
    let hU = (hRho[*, 1..] + hRho[*, ..e1]) / 2.
    let hV = (hRho[1.., *] + hRho[..e2, *]) / 2.

-    let pos: GridPos =
-        {
-            rho = nc[ "lon_rho" ].GetData() :?> float [,], nc[ "lat_rho" ].GetData() :?> float [,]
-            u = nc[ "lon_u" ].GetData() :?> float [,], nc[ "lat_u" ].GetData() :?> float [,]
-            v = nc[ "lon_v" ].GetData() :?> float [,], nc[ "lat_v" ].GetData() :?> float [,]
-        }
-    let angle = nc[ "angle" ].GetData() :?> float [,]
-    let h = { rho = hRho.ToArray(); u = hU.ToArray(); v = hV.ToArray() }
+    let pos: GridPos = {
+        rho = nc["lon_rho"].GetData () :?> float[,], nc["lat_rho"].GetData () :?> float[,]
+        u = nc["lon_u"].GetData () :?> float[,], nc["lat_u"].GetData () :?> float[,]
+        v = nc["lon_v"].GetData () :?> float[,], nc["lat_v"].GetData () :?> float[,]
+    }
+    let angle = nc["angle"].GetData () :?> float[,]
+    let h = { rho = hRho.ToArray (); u = hU.ToArray (); v = hV.ToArray () }
    {
        h = h
-        z = layerZ csR h
+        z = layerZ srho h
        pos = pos
        wetMask = readWetMask nc
        angle = angle
@@ -152,31 +133,28 @@ let readGrid (nc: DataSet) =

 type IRomsToFvcom =
    abstract uv: bool
-    abstract rhoToFvcom: 'a [,] -> 'a []
-    abstract uToFvcom: 'a [,] -> 'a []
-    abstract vToFvcom: 'a [,] -> 'a []
+    abstract rhoToFvcom: 'a[,] -> 'a[]
+    abstract uToFvcom: 'a[,] -> 'a[]
+    abstract vToFvcom: 'a[,] -> 'a[]
    abstract rhoAdjoint: FvcomAdjoint
    abstract uAdjoint: FvcomAdjoint
    abstract vAdjoint: FvcomAdjoint

-type AdjoinedGrid =
-    {
-        h: float []
-        zRho: float [] []
-        zU: float [] []
-        zV: float [] []
-        angle: float []
-        u: float []
-        v: float []
-    }
+type AdjoinedGrid = {
+    h: float[]
+    zRho: float[][]
+    zU: float[][]
+    zV: float[][]
+    angle: float[]
+    u: float[]
+    v: float[]
+}

 let adjoinGirds (r2f: IRomsToFvcom) (roms: RomsGrid) =
    Log.Information "Adjoining grids."
    {
        h = r2f.rhoToFvcom roms.h.rho
-        zRho =
-            roms.z
-            |> Array.map (fun h -> r2f.rhoToFvcom h.rho)
+        zRho = roms.z |> Array.map (fun h -> r2f.rhoToFvcom h.rho)
        zU = roms.z |> Array.map (fun h -> r2f.uToFvcom h.u)
        zV = roms.z |> Array.map (fun h -> r2f.vToFvcom h.v)
        angle = r2f.rhoToFvcom roms.angle
@@ -201,13 +179,13 @@ let genRomsToFvcomAdjoints (proj: IProj) (fvcom: Grid) (roms: RomsGrid) =
        let vM = getNearestCell proj fvcom (roms.pos.v, roms.wetMask.v)
        uM, vM
    { new IRomsToFvcom with
-        member x.uv = true
-        member x.rhoToFvcom m = adjoinMatrix rhoM m
-        member x.uToFvcom m = adjoinMatrix (fst uvMappings) m
-        member x.vToFvcom m = adjoinMatrix (snd uvMappings) m
-        member x.rhoAdjoint = rhoM
-        member x.uAdjoint = fst uvMappings
-        member x.vAdjoint = snd uvMappings
+        member this.uv = true
+        member this.rhoToFvcom m = adjoinMatrix rhoM m
+        member this.uToFvcom m = adjoinMatrix (fst uvMappings) m
+        member this.vToFvcom m = adjoinMatrix (snd uvMappings) m
+        member this.rhoAdjoint = rhoM
+        member this.uAdjoint = fst uvMappings
+        member this.vAdjoint = snd uvMappings
    }

 let adjoinRomsToFvcom (proj: IProj) (fvcom: Grid) (roms: RomsGrid) (props: RomsProps) =
@@ -216,8 +194,7 @@ let adjoinRomsToFvcom (proj: IProj) (fvcom: Grid) (roms: RomsGrid) (props: RomsP

 let computeAngles (grid: FvcomGrid) (adj: AdjoinedGrid) =
    let a = adj.angle
-    grid.Elem
-    |> Array.map (fun (i, j, k) -> (a[i] + a[j] + a[k]) / 3.0)
+    grid.Elem |> Array.map (fun (i, j, k) -> (a[i] + a[j] + a[k]) / 3.0)

 let calcUVBar (fvcom: FvcomGrid) (u, v) =
    let h = fvcom.Bathymetry
@@ -231,23 +208,18 @@ let calcUVBar (fvcom: FvcomGrid) (u, v) =
        [|
            for z = 0 to nZ - 1 do
                fvcom.Elem
-                |> Array.mapi (fun n (i, j, k) ->
-                    hc[n] * (s[z, i] + s[z, j] + s[z, k] |> float)
-                    / 3.0)
+                |> Array.mapi (fun n (i, j, k) -> hc[n] * (s[z, i] + s[z, j] + s[z, k] |> float) / 3.0)
        |]
        |> matrix
    let dz =
-        sigz.ToColumnArrays()
-        |> Array.map (
-            Array.pairwise
-            >> Array.map (fun (x0, x1) -> x1 - x0)
-        )
+        sigz.ToColumnArrays ()
+        |> Array.map (Array.pairwise >> Array.map (fun (x0, x1) -> x1 - x0))
        |> matrix
-        |> fun x -> x.PointwiseAbs().Transpose()
+        |> fun x -> x.PointwiseAbs().Transpose ()
    let bar (m: Matrix<float>) =
        m.PointwiseMultiply dz
        |> fun x ->
-            let s = x.ColumnSums()
+            let s = x.ColumnSums ()
            s.PointwiseDivide hc
    let u' = Matrix.Build.DenseOfRowArrays u
    let v' = Matrix.Build.DenseOfRowArrays v
@@ -256,9 +228,7 @@ let calcUVBar (fvcom: FvcomGrid) (u, v) =
 // Linear interpolation to nearest vertical neighbor
 let mkAllDepthInterpolCoefs (fvcom: Fvcom.FvcomGrid) (roms: AdjoinedGrid) =
    Log.Information "Computing interpolation coefficients."
-    let uv =
-        Array.zip roms.u roms.v
-        |> Array.map (fun (u, v) -> (u + v) * 0.5)
+    let uv = Array.zip roms.u roms.v |> Array.map (fun (u, v) -> (u + v) * 0.5)
    {
        iRho = mkDepthInterpolCoefs fvcom.Siglay roms.h roms.zRho
        iU = mkDepthInterpolCoefs fvcom.SiglayCenter uv roms.zU
--- a/src/Smoothing.fs
+++ b/src/Smoothing.fs
@@ -7,61 +7,61 @@ open ProjNet.FSharp

 let private OOBVAL = -100000f

-let smooth (prop: single []) source =
+let smooth (prop: single[]) source =
    source
    |> Array.fold (fun a x -> a + prop[x]) 0f
    |> fun x -> x / single source.Length

-let smooth' (prop: single []) source =
+let smooth' (prop: single[]) source =
    source
    |> Array.fold (fun (n, a) x -> if prop[x] > OOBVAL then n + 1, a + prop[x] else n, a) (0, 0f)
    |> fun (n, a) -> a / single n

-let smoothNodes (idx: NeighborIndex) (prop: single []) =
-    idx.NodesAroundNode.Values
-    |> Seq.toArray
-    |> Array.Parallel.map (smooth prop)
+let smoothNodes (idx: NeighborIndex) (prop: single[]) =
+    idx.NodesAroundNode.Values |> Seq.toArray |> Array.Parallel.map (smooth prop)

-let smoothNodes2D (idx: NeighborIndex) (prop: single [,]) =
+let smoothNodes2D (idx: NeighborIndex) (prop: single[,]) =
    for k = 0 to Array2D.length1 prop - 1 do
        Log.Debug $"2D smoothing nodes z = {k}"
-        prop[k, *] <- idx.NodesAroundNode.Values
-                      |> Seq.toArray
-                      |> Array.Parallel.map (smooth prop[k, *])
+        prop[k, *] <-
+            idx.NodesAroundNode.Values
+            |> Seq.toArray
+            |> Array.Parallel.map (smooth prop[k, *])
    prop

-let smoothNodes3D (idx: NeighborIndex) (prop: single [,,]) =
+let smoothNodes3D (idx: NeighborIndex) (prop: single[,,]) =
    for k = 0 to Array3D.length1 prop - 1 do
        for n = 0 to Array3D.length2 prop - 1 do
            Log.Debug $"3D smoothing nodes z = {n}"
-            prop[k, n, *] <- idx.NodesAroundNode.Values
-                             |> Seq.toArray
-                             |> Array.Parallel.map (smooth prop[k, n, *])
+            prop[k, n, *] <-
+                idx.NodesAroundNode.Values
+                |> Seq.toArray
+                |> Array.Parallel.map (smooth prop[k, n, *])
    prop

-let smoothElements (idx: NeighborIndex) (grid: Grid) (prop: single []) =
+let smoothElements (idx: NeighborIndex) (grid: Grid) (prop: single[]) =
    grid.Elem
    |> Array.Parallel.mapi (fun n _ ->
        let elems = getElemsSurroundingElem idx grid n
        smooth prop elems)

-let smoothElements2D (idx: NeighborIndex) (grid: Grid) (prop: single [,]) =
+let smoothElements2D (idx: NeighborIndex) (grid: Grid) (prop: single[,]) =
    for k = 0 to Array2D.length1 prop - 1 do
        Log.Debug $"2D smoothing elements z = {k}"
        prop[k, *] <- smoothElements idx grid prop[k, *]
    prop

-let smoothElements3D (idx: NeighborIndex) (grid: Grid) (prop: single [,,]) =
+let smoothElements3D (idx: NeighborIndex) (grid: Grid) (prop: single[,,]) =
    for k = 0 to Array3D.length1 prop - 1 do
        for n = 0 to Array3D.length2 prop - 1 do
            Log.Debug $"3D smoothing elements z = {n}"
            prop[k, n, *] <- smoothElements idx grid prop[k, n, *]
    prop

-let private rectifyOob (getSurrounding: int -> int []) (oob: int []) (prop: single []) =
+let private rectifyOob (getSurrounding: int -> int[]) (oob: int[]) (prop: single[]) =
    for i in oob do
        prop[i] <- OOBVAL // reset
-    let rec rectify (oob: int []) =
+    let rec rectify (oob: int[]) =
        Log.Debug $"rectify oob remaining {oob.Length}"
        for n in oob do
            let ns = getSurrounding n
@@ -73,34 +73,34 @@ let private rectifyOob (getSurrounding: int -> int []) (oob: int []) (prop: sing
    rectify oob
    prop

-let rectifyOutOfBoundsNodes (idx: NeighborIndex) (oob: int []) (prop: single []) =
+let rectifyOutOfBoundsNodes (idx: NeighborIndex) (oob: int[]) (prop: single[]) =
    let f = getNodesSurroundingNode idx
    rectifyOob f oob prop

-let rectifyOutOfBoundsNodes2D (idx: NeighborIndex) (oob: int []) (prop: single [,]) =
+let rectifyOutOfBoundsNodes2D (idx: NeighborIndex) (oob: int[]) (prop: single[,]) =
    let f = getNodesSurroundingNode idx
    for k = 0 to Array2D.length1 prop - 1 do
        prop[k, *] <- rectifyOob f oob prop[k, *]
    prop

-let rectifyOutOfBoundsNodes3D (idx: NeighborIndex) (oob: int []) (prop: single [,,]) =
+let rectifyOutOfBoundsNodes3D (idx: NeighborIndex) (oob: int[]) (prop: single[,,]) =
    let f = getNodesSurroundingNode idx
    for k = 0 to Array3D.length1 prop - 1 do
        for n = 0 to Array3D.length2 prop - 1 do
            prop[k, n, *] <- rectifyOob f oob prop[k, n, *]
    prop

-let rectifyOutOfBoundsElements (idx: NeighborIndex) (grid: Grid) (oob: int []) (prop: single []) =
+let rectifyOutOfBoundsElements (idx: NeighborIndex) (grid: Grid) (oob: int[]) (prop: single[]) =
    let f = getElemsSurroundingElem idx grid
    rectifyOob f oob prop

-let rectifyOutOfBoundsElements2D (idx: NeighborIndex) (grid: Grid) (oob: int []) (prop: single [,]) =
+let rectifyOutOfBoundsElements2D (idx: NeighborIndex) (grid: Grid) (oob: int[]) (prop: single[,]) =
    let f = getElemsSurroundingElem idx grid
    for k = 0 to Array2D.length1 prop - 1 do
        prop[k, *] <- rectifyOob f oob prop[k, *]
    prop

-let rectifyOutOfBoundsElements3D (idx: NeighborIndex) (grid: Grid) (oob: int []) (prop: single [,,]) =
+let rectifyOutOfBoundsElements3D (idx: NeighborIndex) (grid: Grid) (oob: int[]) (prop: single[,,]) =
    let f = getElemsSurroundingElem idx grid
    for k = 0 to Array3D.length1 prop - 1 do
        for n = 0 to Array3D.length2 prop - 1 do
--- a/src/Thredds.fs
+++ b/src/Thredds.fs
@@ -1,6 +1,7 @@
 module Oceanbox.FvcomKit.Thredds

 open System
+
 open Microsoft.Research.Science.Data
 open Microsoft.Research.Science.Data.NetCDF4
 open Serilog
@@ -21,21 +22,20 @@ type Kind =
        | Avg -> "avg"
        | Qck -> "qck"

-let dateRange (start: DateTime) days =
+let dateRange (start: DateTime) days : DateTime list =
    if days < 0 then
-        List.unfold (fun d -> if d < days then None else Some(start.AddDays d, d - 1)) 0
+        List.unfold (fun d -> if d < days then None else Some (start.AddDays d, d - 1)) 0
    else
-        List.unfold (fun d -> if d > days then None else Some(start.AddDays d, d + 1)) 0
+        List.unfold (fun d -> if d > days then None else Some (start.AddDays d, d + 1)) 0

-let tryOpenArchive url =
-    let uri = NetCDFUri()
-    uri.Url <- url
+let tryOpenArchive url : DataSet option =
+    let uri = NetCDFUri ()
+    do uri.Url <- url
    try
        let ds = NetCDFDataSet.Open uri
        Some ds
-    with
-    | e ->
-        Log.Debug e.Message
+    with e ->
+        do Log.Debug e.Message
        None

 let rec tryOpenThredds =
@@ -43,18 +43,22 @@ let rec tryOpenThredds =
    | [] -> None
    | x :: xs ->
        match tryOpenArchive x with
+        | None ->
+            tryOpenThredds xs
        | Some ds ->
-            Log.Debug $"thredds: {x}"
-            Some(x, ds)
-        | None -> tryOpenThredds xs
+            do Log.Debug $"thredds: {x}"
+            Some (x, ds)

 let tryGetTimeIndex (ds: DataSet) (t: DateTime) =
-    let ot = ds[ "ocean_time" ].GetData() :?> double []
-    let t0 = DateTimeOffset.FromUnixTimeSeconds(ot[0] |> int64)
-    let t1 = DateTimeOffset.FromUnixTimeSeconds(ot[^0] |> int64)
-    Log.Debug $"t={t} t0={t0} tn={t1}"
+    let ot = ds["ocean_time"].GetData () :?> double[]
+    let first = Array.head ot
+    let last = Array.last ot
+    let t0 = DateTimeOffset.FromUnixTimeSeconds (first |> int64)
+    let t1 = DateTimeOffset.FromUnixTimeSeconds (last |> int64)
+    do Log.Debug $"t={t} t0={t0} tn={t1}"
+
    if t < t0.DateTime || t > t1.DateTime then
-        Log.Error "time is out of bounds"
+        do Log.Error "time is out of bounds"
        None
    else
        let dt = t - t0.DateTime
--- a/src/Types.fs
+++ b/src/Types.fs
@@ -2,7 +2,26 @@ module Oceanbox.FvcomKit.Types

 open System

-let unzip2D (array: _ [,]) =
+type Vertex = float * float
+
+[<Struct>]
+type BBox = {
+    minX: float
+    maxX: float
+    minY: float
+    maxY: float
+    center: float * float
+} with
+    static member empty = {
+        minX = Double.MaxValue
+        maxX = Double.MinValue
+        minY = Double.MaxValue
+        maxY = Double.MinValue
+        center = 0, 0
+    }
+
+
+let unzip2D (array: _[,]) =
    let len1 = Array2D.length1 array
    let len2 = Array2D.length2 array
    let res1 = Array2D.zeroCreate len1 len2
@@ -14,7 +33,7 @@ let unzip2D (array: _ [,]) =
            res2[i, j] <- y
    res1, res2

-let inline propTo3D (p: 'a [,]) =
+let inline propTo3D (p: 'a[,]) =
    let z = Array2D.length1 p
    let n = Array2D.length2 p
    let m = Array3D.zeroCreate 1 z n
@@ -23,7 +42,7 @@ let inline propTo3D (p: 'a [,]) =
            m[0, i, j] <- p[i, j] |> single
    m

-let inline propTo2D (p: 'a []) =
+let inline propTo2D (p: 'a[]) =
    let n = p.Length
    let m = Array2D.zeroCreate 1 n
    for i = 0 to n - 1 do
--- a/src/default.nix
+++ b/src/default.nix
@@ -0,0 +1,41 @@
+{
+  SDSLite,
+  projnet,
+  dotnet-sdk,
+  dotnet-runtime,
+  buildDotnetModule,
+}:
+let
+  name = "Oceanbox.Fvcomkit";
+  projectFile = ./Oceanbox.FvcomKit.fsproj;
+  versionMatch = builtins.match ".*<Version>([^<]+)</Version>.*" (
+    builtins.readFile projectFile
+  );
+  version = builtins.head versionMatch;
+in
+buildDotnetModule {
+  name = name;
+  pname = name;
+  version = version;
+
+  src = ./.;
+
+  buildInputs = [
+    projnet
+    SDSLite
+  ];
+
+  projectFile = "Oceanbox.FvcomKit.fsproj";
+  inherit
+    dotnet-sdk
+    dotnet-runtime
+    ;
+
+  nugetDeps = ./deps.json;
+
+  packNupkg = true;
+
+  # NOTE(mrtz): Can't package nuget without it
+  # [ref](https://github.com/dotnet/fsharp/issues/12320)
+  dotnetFlags = "--property:TargetsForTfmSpecificContentInPackage=";
+}
--- a/src/deps.json
+++ b/src/deps.json
@@ -0,0 +1,272 @@
+[
+  {
+    "pname": "DynamicInterop",
+    "version": "0.9.1",
+    "hash": "sha256-IB76dA0+K/y/2s/qYL7AfVOF0+6W2hVIBgf9YdZ1oJY="
+  },
+  {
+    "pname": "FSharp.Core",
+    "version": "9.0.303",
+    "hash": "sha256-AxR6wqodeU23KOTgkUfIgbavgbcSuzD4UBP+tiFydgA="
+  },
+  {
+    "pname": "FSharp.Data",
+    "version": "6.3.0",
+    "hash": "sha256-zhVkSfqCljqr6UR0IUMOHUBlR61PvwYKq63PQ09yJPM="
+  },
+  {
+    "pname": "FSharp.Data",
+    "version": "6.4.1",
+    "hash": "sha256-+Z7zbD8cKmhHJWg7Z8XHJ8IeJXCWr/kgRl+VbbsMFw8="
+  },
+  {
+    "pname": "FSharp.Data.Csv.Core",
+    "version": "6.3.0",
+    "hash": "sha256-JdOr3NDmLPohkPpZaWjKqssw0+Wr1lVDtJwTNJ/JhcY="
+  },
+  {
+    "pname": "FSharp.Data.Csv.Core",
+    "version": "6.4.1",
+    "hash": "sha256-oz040beVF7WMONi3n0dPQlZD5deQWnClSXKRijgnw/k="
+  },
+  {
+    "pname": "FSharp.Data.Html.Core",
+    "version": "6.3.0",
+    "hash": "sha256-tSstVvAT9o+0Pr6cIReJOvh0kcthOWgt1CPzgIRoYRQ="
+  },
+  {
+    "pname": "FSharp.Data.Html.Core",
+    "version": "6.4.1",
+    "hash": "sha256-QBbvE8WXUVjS/0mW3aohZBuyfr3M7UGw7kt1oSrlq+s="
+  },
+  {
+    "pname": "FSharp.Data.Http",
+    "version": "6.3.0",
+    "hash": "sha256-/PzzLT0ev4miFswct+YscFDwoaq05BSJATM4fPvxk8o="
+  },
+  {
+    "pname": "FSharp.Data.Http",
+    "version": "6.4.1",
+    "hash": "sha256-0YD/jSCppE1siXrUcTx0OmVdgsjMk+gn0pHp+3GS3V4="
+  },
+  {
+    "pname": "FSharp.Data.Json.Core",
+    "version": "6.3.0",
+    "hash": "sha256-MFe88psxmHWGQYoG8NXi4z33TlWO+dMwOV4NViaUmTM="
+  },
+  {
+    "pname": "FSharp.Data.Json.Core",
+    "version": "6.4.1",
+    "hash": "sha256-0Fmo0f1jC3s+Dime8j2oqLnOK+elqo1xWmktpEYrZlk="
+  },
+  {
+    "pname": "FSharp.Data.Runtime.Utilities",
+    "version": "6.3.0",
+    "hash": "sha256-psc/tsHLYrorjeBBBLviwwA57XMFXUP2ywZqLMzfxac="
+  },
+  {
+    "pname": "FSharp.Data.Runtime.Utilities",
+    "version": "6.4.1",
+    "hash": "sha256-rNo2XQMME1zrPaIezD15P0RoTu8wyhtiJyB99Qp1hcE="
+  },
+  {
+    "pname": "FSharp.Data.WorldBank.Core",
+    "version": "6.3.0",
+    "hash": "sha256-QdL5ylUCvvrhvnnSPWj4MfN7B78hMbb5IRmozK7oJjM="
+  },
+  {
+    "pname": "FSharp.Data.WorldBank.Core",
+    "version": "6.4.1",
+    "hash": "sha256-nUyyziwpY58UnBNpqFoe/1bgDfQIq6gOqtQIwAo7x/c="
+  },
+  {
+    "pname": "FSharp.Data.Xml.Core",
+    "version": "6.3.0",
+    "hash": "sha256-67ftkfQJZ3iD62YKFh8Tu9Fuusb96KlKWxgykP1Wd9U="
+  },
+  {
+    "pname": "FSharp.Data.Xml.Core",
+    "version": "6.4.1",
+    "hash": "sha256-ZiD2aiD5yUZR4CDZl6mh4ji1G3xvm4Yd9Gya/e+D32w="
+  },
+  {
+    "pname": "FSharpPlus",
+    "version": "1.5.0",
+    "hash": "sha256-jQUlF3hsi3xpg+AdTnQw2L+lzbvTh5BIyLXCdVT6u6M="
+  },
+  {
+    "pname": "FSharpPlus",
+    "version": "1.7.0",
+    "hash": "sha256-6hDoDOnMFXQC5Hrk6Fhd+Wj+PbPFXzL9+xLIqgILJuY="
+  },
+  {
+    "pname": "FsPickler",
+    "version": "5.3.2",
+    "hash": "sha256-hjtm55aPJllzcVMPjFP4KYiEEBYtCcrUhbVOR+34agg="
+  },
+  {
+    "pname": "KdTree",
+    "version": "1.4.1",
+    "hash": "sha256-R4+L26pJoliLiwMuxmJDoa3Vf16gBq417fN+iNCy7Yc="
+  },
+  {
+    "pname": "MathNet.Numerics",
+    "version": "5.0.0",
+    "hash": "sha256-RHJCVM6OxquJF7n5Mbe/oNbucBbkge6ULcbAczOgmVo="
+  },
+  {
+    "pname": "MathNet.Numerics.FSharp",
+    "version": "5.0.0",
+    "hash": "sha256-pPbh8JdmMjBgEu84c/qV4YJ+LLr4+c31C6t++u29qBs="
+  },
+  {
+    "pname": "MessagePack",
+    "version": "3.1.3",
+    "hash": "sha256-OBn7iltr/rdE7ZKmv0MCUQSS+6OJKUYtlHdTbhEwzzE="
+  },
+  {
+    "pname": "MessagePack.Annotations",
+    "version": "3.1.3",
+    "hash": "sha256-o+T3u+xaHtW1c7AeWysCmIDUfN8lRhes2LoW5iQBafs="
+  },
+  {
+    "pname": "MessagePackAnalyzer",
+    "version": "3.1.3",
+    "hash": "sha256-5t4Av4CQ8HI7y9aAw+2qcOp+fsY0/3PdaFPJeCEAXQ0="
+  },
+  {
+    "pname": "Microsoft.NET.StringTools",
+    "version": "17.11.4",
+    "hash": "sha256-lWfzY35WQ+iKS9TpuztDTljgF9CIORhFhFEm0p1dVBE="
+  },
+  {
+    "pname": "Microsoft.NETCore.Platforms",
+    "version": "1.1.0",
+    "hash": "sha256-FeM40ktcObQJk4nMYShB61H/E8B7tIKfl9ObJ0IOcCM="
+  },
+  {
+    "pname": "Microsoft.NETCore.Targets",
+    "version": "1.1.0",
+    "hash": "sha256-0AqQ2gMS8iNlYkrD+BxtIg7cXMnr9xZHtKAuN4bjfaQ="
+  },
+  {
+    "pname": "ProjNET",
+    "version": "2.0.0",
+    "hash": "sha256-GjBnuGXmdFagIw9mX51Kpu/nn4gXta6a0cK/dxOWaZY="
+  },
+  {
+    "pname": "runtime.any.System.IO",
+    "version": "4.3.0",
+    "hash": "sha256-vej7ySRhyvM3pYh/ITMdC25ivSd0WLZAaIQbYj/6HVE="
+  },
+  {
+    "pname": "runtime.any.System.Reflection",
+    "version": "4.3.0",
+    "hash": "sha256-ns6f++lSA+bi1xXgmW1JkWFb2NaMD+w+YNTfMvyAiQk="
+  },
+  {
+    "pname": "runtime.any.System.Reflection.Primitives",
+    "version": "4.3.0",
+    "hash": "sha256-LkPXtiDQM3BcdYkAm5uSNOiz3uF4J45qpxn5aBiqNXQ="
+  },
+  {
+    "pname": "runtime.any.System.Runtime",
+    "version": "4.3.0",
+    "hash": "sha256-qwhNXBaJ1DtDkuRacgHwnZmOZ1u9q7N8j0cWOLYOELM="
+  },
+  {
+    "pname": "runtime.any.System.Text.Encoding",
+    "version": "4.3.0",
+    "hash": "sha256-Q18B9q26MkWZx68exUfQT30+0PGmpFlDgaF0TnaIGCs="
+  },
+  {
+    "pname": "runtime.any.System.Threading.Tasks",
+    "version": "4.3.0",
+    "hash": "sha256-agdOM0NXupfHbKAQzQT8XgbI9B8hVEh+a/2vqeHctg4="
+  },
+  {
+    "pname": "runtime.native.System",
+    "version": "4.3.0",
+    "hash": "sha256-ZBZaodnjvLXATWpXXakFgcy6P+gjhshFXmglrL5xD5Y="
+  },
+  {
+    "pname": "runtime.unix.System.Private.Uri",
+    "version": "4.3.0",
+    "hash": "sha256-c5tXWhE/fYbJVl9rXs0uHh3pTsg44YD1dJvyOA0WoMs="
+  },
+  {
+    "pname": "Serilog",
+    "version": "4.2.0",
+    "hash": "sha256-7f3EpCsEbDxXgsuhE430KVI14p7oDUuCtwRpOCqtnbs="
+  },
+  {
+    "pname": "Serilog.Sinks.Console",
+    "version": "6.0.0",
+    "hash": "sha256-QH8ykDkLssJ99Fgl+ZBFBr+RQRl0wRTkeccQuuGLyro="
+  },
+  {
+    "pname": "Serilog.Sinks.File",
+    "version": "6.0.0",
+    "hash": "sha256-KQmlUpG9ovRpNqKhKe6rz3XMLUjkBqjyQhEm2hV5Sow="
+  },
+  {
+    "pname": "Serilog.Sinks.Seq",
+    "version": "9.0.0",
+    "hash": "sha256-NnAkRbxwQGdNXz6DDONRxorNh1nqH2TfAQtokbq5qDw="
+  },
+  {
+    "pname": "System.IO",
+    "version": "4.3.0",
+    "hash": "sha256-ruynQHekFP5wPrDiVyhNiRIXeZ/I9NpjK5pU+HPDiRY="
+  },
+  {
+    "pname": "System.Memory",
+    "version": "4.5.3",
+    "hash": "sha256-Cvl7RbRbRu9qKzeRBWjavUkseT2jhZBUWV1SPipUWFk="
+  },
+  {
+    "pname": "System.Numerics.Vectors",
+    "version": "4.5.0",
+    "hash": "sha256-qdSTIFgf2htPS+YhLGjAGiLN8igCYJnCCo6r78+Q+c8="
+  },
+  {
+    "pname": "System.Private.Uri",
+    "version": "4.3.0",
+    "hash": "sha256-fVfgcoP4AVN1E5wHZbKBIOPYZ/xBeSIdsNF+bdukIRM="
+  },
+  {
+    "pname": "System.Reflection",
+    "version": "4.3.0",
+    "hash": "sha256-NQSZRpZLvtPWDlvmMIdGxcVuyUnw92ZURo0hXsEshXY="
+  },
+  {
+    "pname": "System.Reflection.Emit.ILGeneration",
+    "version": "4.3.0",
+    "hash": "sha256-mKRknEHNls4gkRwrEgi39B+vSaAz/Gt3IALtS98xNnA="
+  },
+  {
+    "pname": "System.Reflection.Emit.Lightweight",
+    "version": "4.3.0",
+    "hash": "sha256-rKx4a9yZKcajloSZHr4CKTVJ6Vjh95ni+zszPxWjh2I="
+  },
+  {
+    "pname": "System.Reflection.Primitives",
+    "version": "4.3.0",
+    "hash": "sha256-5ogwWB4vlQTl3jjk1xjniG2ozbFIjZTL9ug0usZQuBM="
+  },
+  {
+    "pname": "System.Runtime",
+    "version": "4.3.0",
+    "hash": "sha256-51813WXpBIsuA6fUtE5XaRQjcWdQ2/lmEokJt97u0Rg="
+  },
+  {
+    "pname": "System.Text.Encoding",
+    "version": "4.3.0",
+    "hash": "sha256-GctHVGLZAa/rqkBNhsBGnsiWdKyv6VDubYpGkuOkBLg="
+  },
+  {
+    "pname": "System.Threading.Tasks",
+    "version": "4.3.0",
+    "hash": "sha256-Z5rXfJ1EXp3G32IKZGiZ6koMjRu0n8C1NGrwpdIen4w="
+  }
+]
--- a/submodules/FsKdTree
+++ b/submodules/FsKdTree
--- a/submodules/ProjNet.FSharp
+++ b/submodules/ProjNet.FSharp
--- a/test/Tests.fs
+++ b/test/Tests.fs
@@ -1,18 +0,0 @@
-module Tests
-
-open Expecto
-
-let server =
-    testList
-        "Server"
-        [
-            testCase "Adding valid Todo"
-            <| fun _ ->
-                let expectedResult = Ok()
-                Expect.equal (Ok()) expectedResult "Result should be ok"
-        ]
-
-let all = testList "All" [ server ]
-
-[<EntryPoint>]
-let main _ = runTests defaultConfig all
--- a/test/Tests.fsproj
+++ b/test/Tests.fsproj
@@ -1,17 +0,0 @@
-<?xml version="1.0" encoding="utf-8"?>
-<Project Sdk="Microsoft.NET.Sdk">
-  <PropertyGroup>
-    <OutputType>Exe</OutputType>
-    <TargetFramework>net6.0</TargetFramework>
-  </PropertyGroup>
-  <ItemGroup>
-    <Compile Include="Tests.fs" />
-  </ItemGroup>
-  <ItemGroup>
-    <ProjectReference Include="..\src\Oceanbox.FvcomKit.fsproj" />
-  </ItemGroup>
-  <ItemGroup>
-    <PackageReference Include="Expecto" Version="9.0.4" />
-    <PackageReference Update="FSharp.Core" Version="6.0.4" />
-  </ItemGroup>
-</Project>
--- a/xtest/Arome.fs
+++ b/xtest/Arome.fs
@@ -0,0 +1,189 @@
+module Arome
+
+open System
+
+open Microsoft.Research.Science.Data
+
+open Xunit
+open FsUnit.Xunit
+open FsUnit.CustomMatchers
+open Serilog
+open ProjNet.FSharp
+
+open Oceanbox.FvcomKit
+
+let logger =
+    LoggerConfiguration()
+        .MinimumLevel.Is(Events.LogEventLevel.Verbose)
+        .WriteTo.Console()
+        .CreateLogger()
+
+do Log.Logger <- logger
+
+[<Literal>]
+let private path = "/data/archives/Arome/meps_det_sfc_20220102T00Z.nc"
+
+let trans = makeTransform CoordSys.WGS84 (CoordSys.LCCMet ())
+
+// Use Netcdf to open dataset
+let private openDataset (path: string) : DataSet =
+    let sw = Diagnostics.Stopwatch.StartNew()
+    let uri = NetCDF4.NetCDFUri()
+    do uri.FileName <- path
+    do uri.OpenMode <- ResourceOpenMode.ReadOnly
+    do uri.Deflate <- NetCDF4.DeflateLevel.Off
+    let ds = NetCDF4.NetCDFDataSet.Open uri
+    do Log.Debug $"openDataSet: {path} completed in {sw.ElapsedMilliseconds}ms"
+
+    ds
+
+let private aromeGrid = {
+    Arome.SquareGrid.empty with
+        dimensions = 949, 1069
+        BBox = {
+            minX = -18.12241427
+            maxX = 54.24126163
+            minY = 49.7653854
+            maxY = 75.22869642
+            center = 27.12063081, 37.61434821
+        }
+        squareSize = 2.444065489
+}
+
+[<Fact>]
+let ``Create grid``() =
+    use ds = openDataset path
+    let res = Arome.getGrid ds
+
+    //  "Should be able to create square grid from /data/archives/Arome/meps_det_sfc_20220102T00Z.nc"
+    res |> should be (ofCase <@ Result<Arome.SquareGrid, string>.Ok @>)
+
+[<Fact>]
+let ``test grid regularity``() =
+    let trans = makeTransform CoordSys.WGS84 (CoordSys.LCCMet ())
+    use ds = openDataset path
+    let dimensions = ds.Dimensions["x"].Length, ds.Dimensions["y"].Length
+    let longs : float array2d = ds["longitude"].GetData () :?> float[,]
+    let lats : float array2d = ds["latitude"].GetData () :?> float[,]
+    // NOTE: The netcdf file dimensions are defined as (y, x)
+    let width = Array2D.length2 longs
+    let height = Array2D.length1 lats
+    let points : Arome.Points array =
+        let result = Array.create (width * height) Arome.Points.Zero
+        for i in 0 .. height - 1 do
+            for j in 0 .. width - 1 do
+                let lat =  lats[i, j]
+                let lon = longs[i, j]
+                let p = lon, lat
+                // NOTE(simkir): Convert to lambert and undo all I've done... :(
+                let x, y = trans.project p
+                result[i * width + j] <- Arome.Points.OfTuple (single x, single y)
+        result
+
+    let mutable acc = 0.0f
+    for i in 0 .. height - 1 do
+        for j in 0 .. width - 2 do
+            let p0 = points[i * width + j]
+            let p1 = points[i * width + j + 1]
+            let dx = if p0.x < p1.x then p1.x - p0.x else p0.x - p1.x
+
+            acc <- acc + dx
+
+    let count = width * height
+    let avgDiffX = acc / single count
+    Log.Debug("Avg. X distance: {Avg}m", avgDiffX)
+
+    acc <- 0.0f
+    for i in 0 .. height - 2 do
+        for j in 0 .. width - 1 do
+            let p0 = points[i * width + j]
+            let p1 = points[(i + 1) * width + j]
+            let dy = if p0.y < p1.y then p1.y - p0.y else p0.y - p1.y
+
+            acc <- acc + dy
+    let avgDiffY = acc / single count
+    Log.Debug("Avg. Y distance: {Avg}m", avgDiffY)
+
+[<Fact>]
+let ``point within first rect``() =
+    use ds = openDataset path
+    let res = Arome.getGrid ds
+
+    match res with
+    | Ok grid ->
+        let width, _height = grid.dimensions
+        let p0 = grid.points[0 * width + 0]
+        let p1 = grid.points[1 * width + 1]
+        let rect = Arome.Rect.OfPoints p0 p1
+        // { x = -1064984.25f y = -1353489.375f },  { x = -1062474.625f y = -1350985.875f }
+        let p2 = Arome.Points.OfTuple (-1063000.25f, -1353486.0f)
+        let within = Arome.Rect.pointWithin rect p2
+
+        within |> should equal true
+    | Error _ ->
+        failwith "Should create grid"
+
+[<Fact>]
+let ``Read uvs for index (1, 1)``() =
+    use ds = openDataset path
+    let res = Arome.getGrid ds
+    match res with
+    | Ok grid ->
+        let idx =
+            let lat = 50.35
+            let lon =  0.304
+            let p = trans.project((lon, lat))
+            Arome.tryFindIndex grid p
+
+        match idx with
+        | Some (i, j) ->
+            let u, v = Arome.readUV ds 0 i j
+            (u, v) |> should equal (8.530239f, 6.659096f) // "Find uv (8.530239, 6.659096) at index (1, 1)"
+        | None ->
+            failwith "Should find idx"
+    | Error _ ->
+        failwith "Should find grid"
+
+[<Fact>]
+let ``Try search index (1, 1069) based on map coords``() =
+    let lat = 72.800
+    let lon = -17.850
+    let p = trans.project((lon, lat))
+    let res =
+        use ds = openDataset path
+        Arome.getGrid ds
+
+    match res with
+    | Ok grid ->
+        let idx = Arome.tryFindIndex grid p
+
+        idx |> should equal (Some (4, 1067)) // "Should find idx (1, 1069)"
+    | Error _ ->
+        failwith "Should find grid"
+
+[<Fact>]
+let ``Read uvs for pos (60.3562, 5.2178)``() =
+    let lat = 60.3562
+    let lon = 5.2178
+    let p = trans.project((lon, lat))
+    use ds = openDataset path
+    let res = Arome.getGrid ds
+    match res with
+    | Ok grid ->
+        let idx = Arome.tryFindIndex grid p
+        match idx with
+        | Some (i, j) ->
+            let u, v = Arome.readUV ds 0 i j
+            // "Find uv (2.367153168, -0.3955917358) at pos (60.3562, 5.2178) almost in Bergen"
+            (u, v) |> should equal (-0.3750343323f, 7.744056702f)
+        | None ->
+            failwith "Should find idx"
+    | Error _ ->
+        failwith "Should find grid"
+
+[<Fact>]
+let ``Read uvs for t 23 x 43 y 144``() =
+    use ds = openDataset path
+    let u, v = Arome.readUV ds 23 43 144
+    // "Find uv (2.367153168, -0.3955917358) at pos (60.3562, 5.2178) almost in Bergen"
+    (u, v) |> should equal (11.77926254f, 5.764093399f)
--- a/xtest/xtest.fsproj
+++ b/xtest/xtest.fsproj
@@ -0,0 +1,38 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <Nullable>enable</Nullable>
+    <OutputType>Exe</OutputType>
+    <RootNamespace>xtest</RootNamespace>
+    <TargetFramework>net9.0</TargetFramework>
+    <RuntimeIdentifier>linux-x64</RuntimeIdentifier>
+    <!--
+    This template uses native xUnit.net command line options when using 'dotnet run' and
+    VSTest by default when using 'dotnet test'. For more information on how to enable support
+    for Microsoft Testing Platform, please visit:
+      https://xunit.net/docs/getting-started/v3/microsoft-testing-platform
+    -->
+  </PropertyGroup>
+
+  <ItemGroup>
+    <Content Include="xunit.runner.json" CopyToOutputDirectory="PreserveNewest" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Compile Include="Arome.fs" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <PackageReference Include="FsUnit" Version="7.1.1" />
+    <PackageReference Include="FsUnit.xUnit" Version="7.1.1" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="18.0.1" />
+    <PackageReference Include="Oceanbox.SDSLite" Version="2.8.0" />
+    <PackageReference Include="xunit.v3" Version="3.2.1" />
+    <PackageReference Include="xunit.runner.visualstudio" Version="3.1.5" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <ProjectReference Include="..\src\Oceanbox.FvcomKit.fsproj" />
+  </ItemGroup>
+
+</Project>
--- a/xtest/xunit.runner.json
+++ b/xtest/xunit.runner.json
@@ -0,0 +1,3 @@
+{
+    "$schema": "https://xunit.net/schema/current/xunit.runner.schema.json"
+}
@@ -1 +1 @@
 .12.2
 .0.0