Constrain Rdir and Tdir in 2-stream shortwave solver

CliMA · Aug 7, 2023 · ab604ae · ab604ae
1 parent 04aadb1
commit ab604ae
Showing 1 changed file with 9 additions and 7 deletions.
diff --git a/src/rte/RTESolverKernels.jl b/src/rte/RTESolverKernels.jl
@@ -346,6 +346,7 @@ function sw_two_stream!(
     (; τ, ssa, g) = op
     (; Rdif, Tdif, Rdir, Tdir, Tnoscat) = src_sw
     k_min = FT === Float64 ? FT(1e-12) : FT(1e-4)
+    μ₀ = cos(zenith[gcol])
 
     @inbounds for glay in 1:nlay
         # Zdunkowski Practical Improved Flux Method "PIFM"
@@ -368,10 +369,10 @@ function sw_two_stream!(
         @inbounds Tdif[glay, gcol] = RT_term * FT(2) * k * exp_minusktau     # Eqn. 26
 
         # Transmittance of direct, unscattered beam. Also used below
-        @inbounds Tnoscat[glay, gcol] = exp(-τ[glay, gcol] / cos(zenith[gcol]))
+        @inbounds T₀ = Tnoscat[glay, gcol] = exp(-τ[glay, gcol] / cos(zenith[gcol]))
 
         # Direct reflect and transmission
-        k_μ = k * cos(zenith[gcol])
+        k_μ = k * μ₀
         k_γ3 = k * γ3
         k_γ4 = k * γ4
 
@@ -383,23 +384,24 @@ function sw_two_stream!(
             RT_term = ssa[glay, gcol] * RT_term / eps(FT)
         end
 
-        @inbounds Rdir[glay, gcol] =
+        @inbounds Rdir_unconstrained =
             RT_term * (
                 (FT(1) - k_μ) * (α2 + k_γ3) - (FT(1) + k_μ) * (α2 - k_γ3) * exp_minus2ktau -
-                FT(2) * (k_γ3 - α2 * k_μ) * exp_minusktau * Tnoscat[glay, gcol]
+                FT(2) * (k_γ3 - α2 * k_μ) * exp_minusktau * T₀
             )
         #
         # Equation 15, multiplying top and bottom by exp(-k*tau),
         #   multiplying through by exp(-tau/mu0) to
         #   prefer underflow to overflow
         # Omitting direct transmittance
         #
-        @inbounds Tdir[glay, gcol] =
+        @inbounds Tdir_unconstrained =
             -RT_term * (
-                (FT(1) + k_μ) * (α1 + k_γ4) * Tnoscat[glay, gcol] -
-                (FT(1) - k_μ) * (α1 - k_γ4) * exp_minus2ktau * Tnoscat[glay, gcol] -
+                (FT(1) + k_μ) * (α1 + k_γ4) * T₀ - (FT(1) - k_μ) * (α1 - k_γ4) * exp_minus2ktau * T₀ -
                 FT(2) * (k_γ4 + α1 * k_μ) * exp_minusktau
             )
+        Rdir[glay, gcol] = max(FT(0), min(Rdir_unconstrained, (FT(1) - T₀)))
+        Tdir[glay, gcol] = max(FT(0), min(Tdir_unconstrained, (FT(1) - T₀ - Rdir[glay, gcol])))
     end
 end