From 4204c272fb2f9cbccc5b388788404523b33cf992 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= <milankloewer@gmx.de>
Date: Mon, 2 Dec 2024 23:05:16 +0000
Subject: [PATCH] One band longwave radiation implemented (#624)

* One band longwave radiation implemented

* changelog updated

* differential optical thickness corrected

* one-band radiation: indexing bug

* add longwave radiation tests
---
 CHANGELOG.md                      |  2 ++
 src/physics/longwave_radiation.jl | 47 +++++++++++++++++++++++++++++--
 test/longwave_radiation.jl        | 16 +++++++++++
 test/runtests.jl                  |  2 +-
 4 files changed, 64 insertions(+), 3 deletions(-)
 create mode 100644 test/longwave_radiation.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 506864d3f..5aabc0cd8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ## Unreleased
 
+- One-band longwave radiation [#624](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/624)
+- compat entry for FiniteDifferences.jl [#620](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/620)
 - Slab ocean and net surface fluxes [#613](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/613)
 - added compat entries for FiniteDifferences.jl and Enzyme.jl [#620](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/620) [#622](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/622)
 - update to GPUArrays v11, JLArrays v0.2 and remove CUDA from tests [#590](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/590)
diff --git a/src/physics/longwave_radiation.jl b/src/physics/longwave_radiation.jl
index 63021c3cf..a5079d84e 100644
--- a/src/physics/longwave_radiation.jl
+++ b/src/physics/longwave_radiation.jl
@@ -156,12 +156,55 @@ function longwave_radiation!(
     column.outgoing_longwave_radiation = Fₖ
 end
 
-# dummy one band radiation for now
 export OneBandRadiation
 @kwdef struct OneBandRadiation <: AbstractLongwave
     nbands::Int = 1
+    band::Int = 1
 end
 
 OneBandRadiation(SG::SpectralGrid; kwargs...) = OneBandRadiation(; kwargs...)
 initialize!(scheme::OneBandRadiation, model::PrimitiveEquation) = nothing
-longwave_radiation!(column::ColumnVariables, scheme::OneBandRadiation, model::PrimitiveEquation) = nothing
\ No newline at end of file
+
+function longwave_radiation!(
+    column::ColumnVariables{NF},
+    scheme::OneBandRadiation,
+    model::PrimitiveEquation,
+) where NF
+
+    (; nlayers) = column
+    (; band) = scheme                                   # which band is used?
+    τ = view(column.optical_depth_longwave, :, band)    # optical depth of that band
+
+    # precompute differential optical depth dτ
+    dτ = column.a                               # reuse scratch vector a
+    @inbounds for k in eachindex(dτ)
+        dτ[k] = τ[k+1] - τ[k]                   # τ defined on half-levels, dτ is depth across level
+    end
+
+    # precompute Stefan-Boltzmann flux σT^4
+    σ = model.atmosphere.stefan_boltzmann
+    B = column.b                                # reuse scratch vector b
+    (; temp) = column
+    @inbounds for k in eachindex(B, temp)
+        B[k] = σ*temp[k]^4
+    end
+
+    # UPWARD flux U
+    local U::NF = σ*column.surface_temp^4       # boundary condition at surface U(τ=τ(z=0)) = σTₛ⁴
+    column.flux_temp_upward[nlayers+1] += U     # accumulate fluxes
+
+    @inbounds for k in nlayers:-1:1             # integrate from surface up
+        # Frierson et al. 2006, equation 6
+        U -= dτ[k]*(U - B[k])                   # negative because we integrate from surface up in -τ direction
+        column.flux_temp_upward[k] += U         # accumulate that flux
+    end
+
+    # DOWNWARD flux D
+    local D::NF = 0                             # top boundary condition of longwave flux
+                                                # non need to accumulate 0 at top downward flux
+    @inbounds for k in 1:nlayers
+        # Frierson et al. 2006, equation 7
+        D += dτ[k]*(B[k] - D)
+        column.flux_temp_downward[k+1] += D
+    end
+end
\ No newline at end of file
diff --git a/test/longwave_radiation.jl b/test/longwave_radiation.jl
new file mode 100644
index 000000000..348d11117
--- /dev/null
+++ b/test/longwave_radiation.jl
@@ -0,0 +1,16 @@
+@testset "Longwave radiation" begin
+    spectral_grid = SpectralGrid(trunc=31, nlayers=8)
+    @testset for LW in (UniformCooling, JeevanjeeRadiation, OneBandRadiation)
+        longwave_radiation = LW(spectral_grid)
+        model = PrimitiveWetModel(spectral_grid; longwave_radiation)
+        simulation = initialize!(model)
+        run!(simulation, period=Day(5))
+        
+        temp = simulation.diagnostic_variables.grid.temp_grid[:, end]
+
+        for ij in eachindex(temp)
+            # just check that surface temperature isn't completely off
+            @test 200 < temp[ij] < 330
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 6e8bebecd..aa845cbff 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -39,7 +39,7 @@ include("ocean.jl")
 include("large_scale_condensation.jl")
 # include("convection.jl")
 include("optical_depth.jl")
-# include("longwave_radiation.jl")
+include("longwave_radiation.jl")
 # include("shortwave_radiation.jl")
 include("random_process.jl")