lin aerosol

src/CoreRT/CoreKernel/lin_doubling.jl

@@ -62,11 +62,17 @@ function doubling_helper!(pol_type,
         @inbounds @views j₁⁺[:,1,:] .= j₀⁺[:,1,:] .* expk'
         # J⁻₁₂(λ)  = J⁻₀₁(λ).exp(-τ(λ)/μ₀)
         @inbounds @views j₁⁻[:,1,:] .= j₀⁻[:,1,:] .* expk'
-        for ctr=1:4
-            dj₁⁺[ctr,:,1,:] .= (dj₀⁺[ctr,:,1,:] + j₀⁺[:,1,:]*dexpk_fctr) .* expk'
-            # J⁻₁₂(λ)  = J⁻₀₁(λ).exp(-τ(λ)/μ₀)
-            dj₁⁻[ctr,:,1,:] .= (dj₀⁻[ctr,:,1,:] + j₀⁻[:,1,:]*dexpk_fctr) .* expk'
-
+        for ctr=1:3
+            if (ctr==1)
+                dj₁⁺[ctr,:,1,:] .= (dj₀⁺[ctr,:,1,:] + j₀⁺[:,1,:]*dexpk_fctr) .* expk'
+                # J⁻₁₂(λ)  = J⁻₀₁(λ).exp(-τ(λ)/μ₀)
+                dj₁⁻[ctr,:,1,:] .= (dj₀⁻[ctr,:,1,:] + j₀⁻[:,1,:]*dexpk_fctr) .* expk'
+            else
+                dj₁⁺[ctr,:,1,:] .= (dj₀⁺[ctr,:,1,:]) .* expk'
+                # J⁻₁₂(λ)  = J⁻₀₁(λ).exp(-τ(λ)/μ₀)
+                dj₁⁻[ctr,:,1,:] .= (dj₀⁻[ctr,:,1,:]) .* expk'
+            end
+
             # J⁻₀₂(λ) = J⁻₀₁(λ) + T⁻⁻₀₁(λ)[I - R⁻⁺₂₁(λ)R⁺⁻₀₁(λ)]⁻¹[J⁻₁₂(λ) + R⁻⁺₂₁(λ)J⁺₁₀(λ)] (see Eqs.8 in Raman paper draft)
             dj₀⁻[ctr,:,1,:] .= dj₀⁻[ctr,:,1,:] + 
                 (dtt⁺⁺_gp_refl[ctr,:,:,:] ⊠ (j₁⁻ + r⁻⁺ ⊠ j₀⁺)) +
@@ -169,7 +175,7 @@ end
     i = mod(iμ, n_stokes)
     if (i > 2)
         J₀⁻[iμ, 1, n] = - J₀⁻[iμ, 1, n] 
-        dJ₀⁻[1:4, iμ, 1, n] .= - dJ₀⁻[1:4, iμ, 1, n] 
+        dJ₀⁻[1:3, iμ, 1, n] .= - dJ₀⁻[1:3, iμ, 1, n] 
     end
 end
 

src/CoreRT/CoreKernel/lin_elemental.jl

@@ -112,7 +112,7 @@ function elemental!(pol_type, SFI::Bool,
                             dτ::AbstractArray,
                             lin_dτ::AbstractArray,
                             computed_layer_properties::CoreScatteringOpticalProperties,
-                            lin_computed_layer_properties::linCoreScatteringOpticalProperties,
+                            #lin_computed_layer_properties::linCoreScatteringOpticalProperties,
                             m::Int,                     # m: fourier moment
                             ndoubl::Int,                # ndoubl: number of doubling computations needed 
                             scatter::Bool,              # scatter: flag indicating scattering
@@ -183,7 +183,13 @@ function elemental!(pol_type, SFI::Bool,
         dt⁻⁻[1,:,:] .= dt⁺⁺
         dt⁺⁺[2:3,:,:] .= 0
         dt⁻⁻[2:3,:,:] .= 0
-    end    
+    end  
+    dt⁺⁺[1,:,:,:] .*= 2^(-ndoubl)
+    dt⁻⁻[1,:,:,:] .*= 2^(-ndoubl)
+    dr⁻⁺[1,:,:,:] .*= 2^(-ndoubl)
+    dr⁺⁻[1,:,:,:] .*= 2^(-ndoubl)
+    dj₀⁺[1,:,1,:] .*= 2^(-ndoubl)
+    dj₀⁻[1,:,1,:] .*= 2^(-ndoubl)
 end
 
 @kernel function get_elem_rt!(r⁻⁺, t⁺⁺, dr⁻⁺, dt⁺⁺, ϖ_λ, dτ_λ, Z⁻⁺, Z⁺⁺, μ, wct) 
@@ -204,7 +210,7 @@ end
         dr⁻⁺[1,i,j,n] = 
             ϖ_λ[n] * Z⁻⁺[i,j,n2] * 
             (1 / μ[i]) * wct[j] * 
-            tmpM
+            tmpM 
 
         dr⁻⁺[2, i,j,n] = 
             Z⁻⁺[i,j,n2] * 
@@ -263,7 +269,11 @@ end
     nothing
 end
 
-@kernel function get_elem_rt_SFI!(J₀⁺, J₀⁻, dJ₀⁺, dJ₀⁻, ϖ_λ, dτ_λ, τ_sum, Z⁻⁺, Z⁺⁺, μ, ndoubl, wct02, nStokes ,I₀, iμ0, D)
+@kernel function get_elem_rt_SFI!(J₀⁺, J₀⁻, 
+                dJ₀⁺, dJ₀⁻, 
+                ϖ_λ, dτ_λ, τ_sum, Z⁻⁺, Z⁺⁺, 
+                μ, ndoubl, wct02, nStokes,
+                I₀, iμ0, D)
     i_start  = nStokes*(iμ0-1) + 1 
     i_end    = nStokes*iμ0
 
@@ -328,14 +338,14 @@ end
 
     J₀⁺[i, 1, n] *= exp(-τ_sum[n]/μ[i_start]) # how to do this?! Add a fourth derivative to RT kernel elements (only for J terms)
     J₀⁻[i, 1, n] *= exp(-τ_sum[n]/μ[i_start]) # 1: wrt τ, 2: wrt ϖ, 3: wrt Z, 4: wrt τ_sum
-    J₀⁺[4, i, 1, n] = - J₀⁺[i, 1, n]/μ[i_start] 
-    J₀⁻[4, i, 1, n] = - J₀⁻[i, 1, n]/μ[i_start] 
+    #dJ₀⁺[4, i, 1, n] = - J₀⁺[i, 1, n]/μ[i_start] 
+    #dJ₀⁻[4, i, 1, n] = - J₀⁻[i, 1, n]/μ[i_start] 
     if ndoubl >= 1
         J₀⁻[i, 1, n] = D[i,i]*J₀⁻[i, 1, n] #D = Diagonal{1,1,-1,-1,...Nquad times}
         dJ₀⁻[1, i, 1, n] = D[i,i]*dJ₀⁻[1, i, 1, n]
         dJ₀⁻[2, i, 1, n] = D[i,i]*dJ₀⁻[2, i, 1, n]
         dJ₀⁻[3, i, 1, n] = D[i,i]*dJ₀⁻[3, i, 1, n]
-        dJ₀⁻[4, i, 1, n] = D[i,i]*dJ₀⁻[4, i, 1, n]
+        #dJ₀⁻[4, i, 1, n] = D[i,i]*dJ₀⁻[4, i, 1, n]
     end  
     nothing
 end
@@ -374,7 +384,7 @@ end
     if ndoubl>1
         if mod(i, pol_n) > 2
             J₀⁻[i, 1, n] = - J₀⁻[i, 1, n]
-            dJ₀⁻[1:4, i, 1, n] .= - dJ₀⁻[1:4, i, 1, n]
+            dJ₀⁻[1:3, i, 1, n] .= - dJ₀⁻[1:3, i, 1, n]
         end 
     end
     nothing

src/CoreRT/CoreKernel/lin_rt_kernel.jl

@@ -14,8 +14,8 @@ function rt_kernel!(RS_type::noRS{FT},
                     computed_layer_properties::M, 
                     lin_computed_layer_properties, 
                     scattering_interface, 
-                    τ_sum, 
-                    lin_τ_sum,
+                    τ_sum, μ₀,
+                    #lin_τ_sum,
                     m, quad_points, 
                     I_static, 
                     architecture, 
@@ -105,16 +105,25 @@ function rt_kernel!(RS_type::noRS{FT},
 
             lin_composite_layer.dJ₀⁺[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁺[1,:,1,:] +
                                                 lin_ϖ[ctr]*lin_added_layer.dj₀⁺[2,:,1,:] +
-                                                lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dj₀⁺[3,:,1,:] +
-                                                lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁺[4,:,1,:]
+                                                lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dj₀⁺[3,:,1,:] #+
+                                                #lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁺[4,:,1,:]
 
             lin_composite_layer.dJ₀⁻[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁻[1,:,1,:] +
                                                 lin_ϖ[ctr]*lin_added_layer.dj₀⁻[2,:,1,:] +
-                                                lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dj₀⁻[3,:,1,:] +
-                                                lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁻[4,:,1,:]
+                                                lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dj₀⁻[3,:,1,:] #+
+                                                #lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁻[4,:,1,:]
         end
+        otmpE = exp.(-τ_sum/μ₀)
+        oτ_sum = τ_sum
+        olin_τ = lin_τ
+        odtmpE = 0.0 .* lin_τ
+        olin_τ = lin_τ
         # If this is not the TOA, perform the interaction step
     else        
+        tmpE = exp.(-τ_sum/μ₀)
+        composite_layer.J₀⁺[:], composite_layer.J₀⁻[:] = (added_layer.j₀⁺.*tmpE, added_layer.j₀⁻.*tmpE)
+
+
         for ctr=1:nparams
             lin_added_layer.dxt⁺⁺[ctr,:,:,:] = lin_τ[ctr,:].*lin_added_layer.dt⁺⁺[1,:,:,:] +
                                                 lin_ϖ[ctr]*lin_added_layer.dt⁺⁺[2,:,:,:] + 
@@ -132,21 +141,31 @@ function rt_kernel!(RS_type::noRS{FT},
                                                 lin_ϖ[ctr]*lin_added_layer.dr⁺⁻[2,:,:,:] +
                                                 lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dr⁺⁻[3,:,:,:]
 
+
+            dtmpE[ctr,:] = exp.(-(τ_sum-oτ_sum)/μ₀) .* (odtmpE[ctr,:] .- otmpE .* olin_τ[ctr,:]./μ₀)
+
             lin_added_layer.dxj₀⁺[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁺[1,:,1,:] +
                                                 lin_ϖ[ctr]*lin_added_layer.dj₀⁺[2,:,1,:] +
-                                                lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dj₀⁺[3,:,1,:] +
-                                                lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁺[4,:,1,:]
-
+                                                lin_Z⁺⁺[ctr,:,:,:].*lin_added_layer.dj₀⁺[3,:,1,:] #+
+                                                #lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁺[4,:,1,:]
+            lin_added_layer.dxj₀⁺[ctr,:,1,:] = lin_added_layer.dxj₀⁺[ctr,:,1,:] .* tmpE .+
+                                                added_layer.j₀⁺.* dtmpE[ctr,:]
             lin_added_layer.dxj₀⁻[ctr,:,1,:] = lin_τ[ctr,:].*lin_added_layer.dj₀⁻[1,:,1,:] +
                                                 lin_ϖ[ctr]*lin_added_layer.dj₀⁻[2,:,1,:] +
-                                                lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dj₀⁻[3,:,1,:] +
-                                                lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁻[4,:,1,:]
+                                                lin_Z⁻⁺[ctr,:,:,:].*lin_added_layer.dj₀⁻[3,:,1,:] #+
+                                                #lin_τ_sum[ctr,:,:,:].*lin_added_layer.dj₀⁻[4,:,1,:]
+            lin_added_layer.dxj₀⁻[ctr,:,1,:] = lin_added_layer.dxj₀⁻[ctr,:,1,:] .* tmpE .+
+                                                added_layer.j₀⁻ .* dtmpE[ctr,:]
         end
         @timeit "interaction" interaction!(RS_type, scattering_interface, SFI, composite_layer, lin_composite_layer, added_layer, lin_added_layer, I_static)
+        otmpE = tmpE
+        oτ_sum = τ_sum
+        odtmpE = dtmpE
+        olin_τ = lin_τ
     end
 end
 
-#=
+
 function get_dtau_ndoubl(computed_layer_properties::CoreScatteringOpticalProperties,
         lin_computed_layer_properties::linCoreScatteringOpticalProperties, 
         quad_points::QuadPoints{FT}) where {FT}
@@ -157,7 +176,7 @@ function get_dtau_ndoubl(computed_layer_properties::CoreScatteringOpticalPropert
     _, ndoubl = doubling_number(dτ_max, maximum(τ .* ϖ))
     # Compute dτ vector
     dτ = τ ./ 2^ndoubl
-    lin_dτ = lin_τ / 2^ndoubl
+    lin_dτ = lin_τ ./ 2^ndoubl
 
     return dτ, lin_dτ, ndoubl
 end
@@ -173,7 +192,7 @@ function get_dtau_ndoubl(computed_layer_properties::CoreDirectionalScatteringOpt
     _, ndoubl = doubling_number(dτ_max, maximum(τ .* ϖ))
     # Compute dτ vector
     dτ = τ ./ 2^ndoubl
-    lin_dτ = lin_τ / 2^ndoubl
+    lin_dτ = lin_τ ./ 2^ndoubl
 
     return dτ, lin_dτ, ndoubl
 end
@@ -191,9 +210,9 @@ function init_layer(computed_layer_properties::CoreDirectionalScatteringOpticalP
     gfct = Array(G)[iμ₀]
     expk = exp.(-dτ*gfct/μ₀)
     #lin_expk = -expk.*lin_dτ*gfct/μ₀
-    dexpk_fctr = -gfct/μ₀
+    lin_expk_fctr = -gfct/μ₀ #.+zeros(length(expk))
 
-    return dτ, lin_dτ, ndoubl, arr_type(expk), FT(dexpk_fctr) #arr_type(lin_expk)
+    return dτ, lin_dτ, ndoubl, arr_type(expk), FT(lin_expk_fctr)#, arr_type(lin_expk)
 end
 
 function init_layer(computed_layer_properties::CoreScatteringOpticalProperties, 
@@ -207,12 +226,12 @@ function init_layer(computed_layer_properties::CoreScatteringOpticalProperties,
                             quad_points)
     expk = exp.(-dτ/μ₀)
     #lin_expk = -expk.*lin_dτ/μ₀
-    dexpk_fctr = -1/μ₀
+    lin_expk_fctr = -1/μ₀
 
-    return dτ, lin_dτ, ndoubl, arr_type(expk), FT(dexpk_fctr)#arr_type(lin_expk)
+    return dτ, lin_dτ, ndoubl, arr_type(expk), FT(lin_expk_fctr)#arr_type(lin_expk)
 end
 
-=#
+#=
 function rt_kernel!(RS_type::Union{RRS{FT}, VS_0to1{FT}, VS_1to0{FT}}, 
         pol_type, SFI, 
         added_layer, lin_added_layer,
@@ -409,5 +428,5 @@ function rt_kernel!(
                 added_layer, lin_added_layer, I_static)
     end
 end
-
+=#
 

src/CoreRT/CoreRT.jl

@@ -87,6 +87,7 @@ include("gpu_batched.jl")                   # Batched operations
 
 # Utilities / Helper Functions
 include("atmo_prof.jl")                     # Helper Functions for Handling Atmospheric Profiles
+include("lin_atmo_prof.jl")  
 include("rt_helper_functions.jl")           # Miscellaneous Utility Functions
 include("rt_set_streams.jl")                # Set streams before RT
 include("parameters_from_yaml.jl")          # Loading in parameters from YAML file

src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl

@@ -132,7 +132,7 @@ function getG_atSun(lod::CoreDirectionalScatteringOpticalProperties,quad_points:
 end
 
 
-function (in::CoreScatteringOpticalProperties, arr_type)
+function expandOpticalProperties(in::CoreScatteringOpticalProperties, arr_type)
     @unpack τ, ϖ, Z⁺⁺, Z⁻⁺ = in 
     @assert length(τ) == length(ϖ) "τ and ϖ sizes need to match"
     if size(Z⁺⁺,3) == 1