From 7dadb87604057b4d8bc09e835e04840d0ce80e7d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= <milankloewer@gmx.de>
Date: Mon, 6 Jan 2025 21:33:16 -0100
Subject: [PATCH] Clamp in ZonalWind to [0, 1] for sqrt (#649)

* clamp in ZonalWind to [0, 1] for sqrt

* update changelog

* Test initialize + time step for uncommon resolutions
---
 CHANGELOG.md                       | 1 +
 src/dynamics/initial_conditions.jl | 2 +-
 test/resolutions.jl                | 9 +++++++++
 test/runtests.jl                   | 1 +
 4 files changed, 12 insertions(+), 1 deletion(-)
 create mode 100644 test/resolutions.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 74bc44421..746ef97cb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+- Bug: ZonalWind had sqrt of negative for uncommon resolutions [#649](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/649)
 - Spectral Gradients are now differentiable, with correctness check in extended CI [#638](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/638)
 - Tracer advection [#579](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/579)
 - Buildkite CI with dummy pipeline [#646](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/646)
diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl
index e7133d43f..0891c6b8d 100644
--- a/src/dynamics/initial_conditions.jl
+++ b/src/dynamics/initial_conditions.jl
@@ -289,7 +289,7 @@ function initialize!(   progn::PrognosticVariables{NF},
             R = radius*perturb_radius # spatial extent of perturbation
 
             # great circle distance to perturbation
-            X = sinφc*sinφ + cosφc*cosφ*cosd(λij-λc)
+            X = clamp(sinφc*sinφ + cosφc*cosφ*cosd(λij-λc), 0, 1)
             X_norm = 1/sqrt(1-X^2)
             r = radius*acos(X)
             exp_decay = exp(-(r/R)^2)
diff --git a/test/resolutions.jl b/test/resolutions.jl
new file mode 100644
index 000000000..0956839cf
--- /dev/null
+++ b/test/resolutions.jl
@@ -0,0 +1,9 @@
+@testset "Single time step at uncommon resolutions" begin
+    @testset for trunc in rand(15:200, 10)
+        spectral_grid = SpectralGrid(; trunc, Grid=FullClenshawGrid)
+        model = PrimitiveWetModel(spectral_grid)
+        simulation = initialize!(model)
+        run!(simulation, period=Day(0))
+        @test simulation.model.feedback.nars_detected == false
+    end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 76f6e7925..b601e54f2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -52,6 +52,7 @@ include("random_process.jl")
 include("stochastic_physics.jl")
 
 # INITIALIZATION AND INTEGRATION
+include("resolutions.jl")
 include("run_speedy.jl")
 
 # EXTENSION