From 58e1ce11c33b95f90801cabeed6a500a945c4327 Mon Sep 17 00:00:00 2001 From: sunitisanghavi Date: Sun, 25 Feb 2024 14:00:32 -0800 Subject: [PATCH] lin aerosol --- src/CoreRT/CoreKernel/lin_doubling.jl | 18 +- src/CoreRT/CoreKernel/lin_elemental.jl | 26 +- src/CoreRT/CoreKernel/lin_rt_kernel.jl | 59 ++- src/CoreRT/CoreRT.jl | 1 + .../compEffectiveLayerProperties.jl | 2 +- .../lin_compEffectiveLayerProperties.jl | 248 +++++++-- src/CoreRT/atmo_prof.jl | 43 +- src/CoreRT/lin_atmo_prof.jl | 486 ++++++++++++++++++ src/CoreRT/lin_model_from_parameters.jl | 125 ++++- src/CoreRT/lin_types.jl | 247 +++++---- src/CoreRT/parameters_from_yaml.jl | 14 +- src/CoreRT/rt_helper_functions.jl | 217 ++++---- src/CoreRT/rt_run.jl | 36 +- src/Scattering/Scattering.jl | 1 + src/Scattering/lin_compute_NAI2.jl | 34 +- src/Scattering/lin_mie_helper_functions.jl | 21 +- src/Scattering/lin_types.jl | 2 +- src/vSmartMOM.jl | 2 +- test/prototyping/Balsamic_OCO2_test.jl | 69 ++- test/prototyping/runner.jl | 88 +++- test/test_parameters/3BandParameters.yaml | 25 +- 21 files changed, 1353 insertions(+), 411 deletions(-) create mode 100644 src/CoreRT/lin_atmo_prof.jl diff --git a/src/CoreRT/CoreKernel/lin_doubling.jl b/src/CoreRT/CoreKernel/lin_doubling.jl index 6ad510ea..4cfb2dff 100644 --- a/src/CoreRT/CoreKernel/lin_doubling.jl +++ b/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 diff --git a/src/CoreRT/CoreKernel/lin_elemental.jl b/src/CoreRT/CoreKernel/lin_elemental.jl index 357fc2c8..758742bf 100644 --- a/src/CoreRT/CoreKernel/lin_elemental.jl +++ b/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 diff --git a/src/CoreRT/CoreKernel/lin_rt_kernel.jl b/src/CoreRT/CoreKernel/lin_rt_kernel.jl index a6f03243..e97918e8 100644 --- a/src/CoreRT/CoreKernel/lin_rt_kernel.jl +++ b/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 - +=# diff --git a/src/CoreRT/CoreRT.jl b/src/CoreRT/CoreRT.jl index a4ba8116..d4d32258 100644 --- a/src/CoreRT/CoreRT.jl +++ b/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 diff --git a/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl b/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl index 6d0d2c28..6681d86d 100644 --- a/src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl +++ b/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 diff --git a/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl b/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl index 1796c430..588d3a60 100644 --- a/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl +++ b/src/CoreRT/LayerOpticalProperties/lin_compEffectiveLayerProperties.jl @@ -1,8 +1,12 @@ +# Construct τ, ϖ, Z⁺⁺, Z⁻⁺ and +# lin_τ, lin_ϖ, linZ⁺⁺, linZ⁻⁺ wrt all state vector parameters +# for each atmospheric layer + function constructCoreOpticalProperties( RS_type::AbstractRamanType{FT}, iBand, m, model, lin) where FT @unpack τ_rayl, τ_aer, τ_abs, aerosol_optics, greek_rayleigh = model - @unpack lin_τ_aer, lin_τ_abs, lin_aerosol_optics = lin + @unpack lin_τ_abs, lin_τ_rayl, lin_τ_aer_psurf, lin_τ_aer_z₀, lin_τ_aer_σz, lin_aerosol_optics = lin #@show typeof(τ_rayl[1]), typeof(τ_aer[1]), typeof(τ_abs[1]) @assert all(iBand .≤ length(τ_rayl)) "iBand exceeded number of bands" @@ -10,13 +14,22 @@ function constructCoreOpticalProperties( arr_type = array_type(model.params.architecture) pol_type = model.params.polarization_type - + #nparams = 7*Naer + 6 + 2 + 1 + 2*n_bands + 2 + # 7 for each aerosol type , + # 6 for CO2 VMR (VMR(z)=a₁f₁(z)+a₂f₂(z)+a₃f₃(z)), + # 2 for H2O VMR (VMR(z)=b₁f₁(z)+b₂f₂(z)), + # 1 for p_surf, + # nbands for ρ₀ (RPV), + # nbands for ρc (RPV), + # 1 for Θ (RPV), + # 1 for k (RPV) + nparams = 23 #excluding surface for now # Quadrature points: μ = Array(model.quad_points.qp_μ ) N = length(model.quad_points.qp_μN) # Number of Aerosols: nAero = size(τ_aer[iBand[1]],1) - nZ = size(τ_rayl[1],2) + Nz = size(τ_rayl[1],2) # Rayleigh Z matrix: Rayl𝐙⁺⁺, Rayl𝐙⁻⁺ = Scattering.compute_Z_moments(pol_type, μ, greek_rayleigh, m, @@ -28,12 +41,49 @@ function constructCoreOpticalProperties( for iB in iBand # Define G here as ones: #G = arr_type(ones(FT,N) - rayl = [CoreScatteringOpticalProperties( - arr_type(τ_rayl[iB][:,i]),RS_type.ϖ_Cabannes[iB], - Rayl𝐙⁺⁺, Rayl𝐙⁻⁺) for i=1:nZ] + #rayl = [CoreScatteringOpticalProperties( + # arr_type(τ_rayl[iB][:,i]),RS_type.ϖ_Cabannes[iB], + # Rayl𝐙⁺⁺, Rayl𝐙⁻⁺) for i=1:Nz] + + #@show size(τ_rayl[iB][:,end]), size(rayl[end].τ) + #@show RS_type.ϖ_Cabannes[iB], rayl[end].ϖ + #@show size(Rayl𝐙⁺⁺), size(rayl[end].Z⁺⁺) + dtmp_τ = [zeros(FT, nparams, size(lin_τ_rayl[1],1)) for i=1:Nz] + for iz = 1:Nz + @show size(lin_τ_rayl[iB],1) + dtmp_τ[iz][1,:] .= lin_τ_rayl[iB][:,iz] + end + dtmp_ϖ = zeros(FT, nparams) + dtmp_Z⁺⁺ = zeros(FT, nparams, size(Rayl𝐙⁺⁺,1), size(Rayl𝐙⁺⁺,1)) + dtmp_Z⁻⁺ = zeros(FT, nparams, size(Rayl𝐙⁺⁺,1), size(Rayl𝐙⁺⁺,1)) + + + #lin_rayl = [linCoreScatteringOpticalProperties( + # arr_type(dtmp_τ[iz]), + # dtmp_ϖ, + # dtmp_Z⁺⁺, dtmp_Z⁻⁺) + # for iz=1:Nz] + rayl = [scatt_paired_xdx(CoreScatteringOpticalProperties( + arr_type(τ_rayl[iB][:,i]),RS_type.ϖ_Cabannes[iB], + Rayl𝐙⁺⁺, Rayl𝐙⁻⁺), + linCoreScatteringOpticalProperties( + arr_type(dtmp_τ[i]), + dtmp_ϖ, + dtmp_Z⁺⁺, dtmp_Z⁻⁺)) + for i=1:Nz] + @show size(rayl) + + #@show size(dtmp_τ[end]), size(lin_rayl[end].lin_τ) + #@show dtmp_ϖ, lin_rayl[end].lin_ϖ + #@show size(dtmp_Z⁺⁺), size(lin_rayl[end].lin_Z⁺⁺) # Initiate combined properties with rayleigh - combo = rayl - lin_combo = [] + combo = rayl #scatt_paired_xdx(rayl, lin_rayl) + @show size(combo[end].x.τ),size(combo[end].dx.lin_τ) + #combo, lin_combo = rayl, lin_rayl + + #@show size(combo.x[1].τ) + #@show size(combo.dx[1].lin_τ) + #lin_combo = [] # Loop over all aerosol types: for i=1:nAero @@ -46,25 +96,48 @@ function constructCoreOpticalProperties( dAerZ⁺⁺ = zeros(4, size(AerZ⁺⁺,1), size(AerZ⁺⁺,2)) dAerZ⁻⁺ = zeros(4, size(AerZ⁺⁺,1), size(AerZ⁺⁺,2)) - for ctr=1:4 + for ctr=1:4 dAerZ⁺⁺[ctr,:,:], dAerZ⁻⁺[ctr,:,:] = Scattering.compute_Z_moments( pol_type, μ, lin_aerosol_optics[iB][i].d_greek_coefs[ctr], m, arr_type=arr_type) end + #@show typeof(AerZ⁺⁺), typeof(aerosol_optics[iB][i]), typeof(FT.(τ_aer[iB][i,:])) # Generate Core optical properties for Aerosols i - aer, lin_aer = createAero.(τ_aer[iB][i,:], - [aerosol_optics[iB][i]], - [AerZ⁺⁺], [AerZ⁻⁺], lin_aerosol_optics[iB][i], dAerZ⁺⁺, dAerZ⁻⁺) - #@show typeof(aer), typeof(combo) + #aer, lin_aer = + aer = createAero(τ_aer[iB][i,:], + (aerosol_optics[iB][i]), + AerZ⁺⁺, AerZ⁻⁺, + nparams, Nz, i, + lin_τ_aer_psurf[iB][i,:], + lin_τ_aer_z₀[iB][i,:], + lin_τ_aer_σz[iB][i,:], + lin_aerosol_optics[iB][i], + dAerZ⁺⁺, dAerZ⁻⁺) + + #combo_aer = scatt_paired_xdx(aer, lin_aer) + #@show size(τ_aer[iB][i,end]) + #@show size(aer), size(aer[end].τ) + #@show aerosol_optics[iB][i].ω̃, aer[end].ϖ + #@show size(AerZ⁺⁺), size(aer[end].Z⁺⁺) + + #@show size(dtmp_τ[end]), size(lin_aer[end].lin_τ) + #@show lin_aerosol_optics[iB][i].dω̃, lin_aer[end].lin_ϖ + #@show size(dAerZ⁺⁺), size(lin_aer[end].lin_Z⁺⁺) + # aer -> CoreScatteringOpticalProperties(τ_mod, ϖ_mod,AerZ⁺⁺, AerZ⁻⁺), + # lin_aer -> linCoreScatteringOpticalProperties(lin_τ_mod, lin_ϖ_mod, dAerZ⁺⁺, dAerZ⁻⁺) + + #@show typeof(aer), typeof(combo) # Mix with previous Core Optical Properties - (combo, lin_combo) = (combo, lin_combo) .+ (aer, lin_aer) + + #lin_combo = lin_combo .+ lin_aer + combo = combo .+ aer + #combo = combo + combo_aer + #(combo, lin_combo) = (combo, lin_combo) .+ (aer, lin_aer) end #@show typeof(combo) - # TODO Type check τ_abs, τ_aer, rayl[i].τ ./ combo[i].τ - # Somewhere here we can add canopy later as well! ### # fScattRayleigh: @@ -72,7 +145,7 @@ function constructCoreOpticalProperties( #@show size(combo) #fScattRayleigh = [FT.(Array(rayl[i].τ ./ combo[i].τ)) for i=1:nZ] fScattRayleigh = [Array(rayl[i].τ ./ combo[i].τ) for i=1:nZ] - lin_fScattRayleigh = [Array(-rayl[i].τ .* lin_combo[i].lin_τ ./ combo[i].τ^2) for i=1:nZ] + #lin_fScattRayleigh = [Array(-rayl[i].τ .* lin_combo[i].lin_τ ./ combo[i].τ^2) for i=1:nZ] #@show fScattRayleigh, rayl[1].τ, combo[1].τ # Create Core Optical Properties merged with trace gas absorptions: @@ -80,15 +153,26 @@ function constructCoreOpticalProperties( #[CoreAbsorptionOpticalProperties(arr_type((τ_abs[iB][:,i]))) for i=1:nZ]) # Gaseous Absorption - gabs = [CoreAbsorptionOpticalProperties( - arr_type(τ_abs[iB][:,i])) for i=1:nZ] - lin_gabs = [linCoreAbsorptionOpticalProperties( - arr_type(lin,τ_abs[iB][:,:,i])) for i=1:nZ] - (combo, lin_combo) = (combo, lin_combo) .+ (gabs, lin_gabs) - push!(band_layer_props, combo) - push!(lin_band_layer_props, lin_combo) + gabs = [abs_paired_xdx(CoreAbsorptionOpticalProperties( + arr_type(τ_abs[iB][:,i])), + linCoreAbsorptionOpticalProperties( + arr_type(lin_τ_abs[iB][:,:,i]))) + for i=1:nZ] + #lin_gabs = [linCoreAbsorptionOpticalProperties( + # arr_type(lin_τ_abs[iB][:,:,i])) for i=1:nZ] + #combo_gabs = abs_paired_xdx(gabs, lin_gabs) + #combo = combo .+ combo_gabs + combo = combo .+ gabs + #(combo, lin_combo) = (combo, lin_combo) .+ (gabs, lin_gabs) + + # Initiate combined properties with rayleigh + #(combo, lin_combo) = (combo, lin_combo) .+ (gabs, lin_gabs) + + push!(band_layer_props, combo.x) + push!(lin_band_layer_props, combo.dx) + #push!(lin_band_layer_props, lin_combo) push!(band_fScattRayleigh,fScattRayleigh) - push!(lin_band_fScattRayleigh,lin_fScattRayleigh) + #push!(lin_band_fScattRayleigh,lin_fScattRayleigh) end layer_opt = [] @@ -98,29 +182,88 @@ function constructCoreOpticalProperties( for iz = 1:nZ push!(layer_opt, prod([band_layer_props[i][iz] for i=1:length(iBand)])); push!(fscat_opt, [band_fScattRayleigh[i][iz] for i=1:length(iBand)]); - push!(lin_layer_opt, prod([lin_band_layer_props[i][:,iz] for i=1:length(iBand)])); - push!(lin_fscat_opt, [lin_band_fScattRayleigh[i][:,iz] for i=1:length(iBand)]); + push!(lin_layer_opt, prod([lin_band_layer_props[i][iz] for i=1:length(iBand)])); + #push!(lin_fscat_opt, [lin_band_fScattRayleigh[i][:,iz] for i=1:length(iBand)]); end # For now just one band_fScattRayleigh - return layer_opt, fscat_opt, lin_layer_opt, lin_fscat_opt + return layer_opt, fscat_opt, lin_layer_opt#, lin_fscat_opt end -function createAero(τAer, aerosol_optics, AerZ⁺⁺, AerZ⁻⁺, lin_aerosol_optics, dAerZ⁺⁺, dAerZ⁻⁺) +function createAero(τAer, + aerosol_optics, + AerZ⁺⁺, AerZ⁻⁺, + nparams, Nz, iAer, + lin_τ_aer_psurf, + lin_τ_aer_z₀, + lin_τ_aer_σz, + lin_aerosol_optics, + dAerZ⁺⁺, dAerZ⁻⁺) + @unpack k_ref, k, fᵗ, ω̃ = aerosol_optics @unpack dk_ref, dk, dfᵗ, dω̃ = lin_aerosol_optics τ_mod = (1-fᵗ * ω̃ ) * τAer; ϖ_mod = (1-fᵗ) * ω̃/(1-fᵗ * ω̃) - - lin_τ_mod[1] = (τ_mod/τAer)*(k/k_ref) + FT = eltype(ϖ_mod) + lin_τ_mod = zeros(5,length(τ_mod)) + lin_ϖ_mod = zeros(5) + lin_τ_mod[1,:] = (τ_mod./τAer)*(k/k_ref) lin_ϖ_mod[1] = 0 + lin_τ_aer_psurf .*= (τ_mod./τAer) + lin_τ_aer_z₀ .*= (τ_mod./τAer) + lin_τ_aer_σz .*= (τ_mod./τAer) for ctr=2:5 #ctr=1 corresponds to the derivative with respect to τ_ref, the rest for the microphysical parameters nᵣ, nᵢ, r₀. and σ₀ mctr = ctr-1 - lin_τ_mod[ctr] = (τ_mod/k)*(dk[mctr] - dk_ref[mctr]*(k/k_ref)) - - τAer*(fᵗ*dω̃[mctr] + dfᵗ[mctr]*ω̃) - lin_ϖ_mod[ctr] = (dω̃[mctr]*(1-fᵗ) - dfᵗ[mctr]*ω̃*(1-ω̃))/(1-fᵗ * ω̃)^2 + lin_τ_mod[ctr,:] = τ_mod * + (dk[mctr]/k - dk_ref[mctr]/k_ref) + - τAer * (fᵗ*dω̃[mctr] + dfᵗ[mctr]*ω̃) + lin_ϖ_mod[ctr] = (dω̃[mctr] - + (dfᵗ[mctr]*ω̃ + fᵗ*dω̃[mctr])*(1-ϖ_mod))/ + (1-fᵗ * ω̃) + end + dtmp_τ = [zeros(FT, nparams) for i=1:Nz] + dtmp_ϖ = zeros(FT, nparams) + dtmp_Z⁺⁺ = zeros(FT, nparams, size(AerZ⁺⁺,1), size(AerZ⁺⁺,1)) + dtmp_Z⁻⁺ = zeros(FT, nparams, size(AerZ⁺⁺,1), size(AerZ⁺⁺,1)) + @show size(τ_mod[end]), length(τ_mod), size(ϖ_mod), size(lin_τ_mod), size(lin_ϖ_mod) + @show nparams, Nz + for iz = 1:Nz + dtmp_τ[iz][1,:] .= lin_τ_aer_psurf[iz] + ctr = 9 + (iAer-1)*7 + dtmp_τ[iz][ctr+6,:] .= lin_τ_aer_z₀[iz] + dtmp_τ[iz][ctr+7,:] .= lin_τ_aer_σz[iz] + dtmp_τ[iz][ctr+1,:] .= lin_τ_mod[1,iz] + dtmp_ϖ[ctr+1] = lin_ϖ_mod[1] + dtmp_Z⁺⁺[ctr+1,:,:] .= 0 + dtmp_Z⁻⁺[ctr+1,:,:] .= 0 + for i=2:5 + mctr=i-1 + dtmp_τ[iz][ctr+i,:] .= lin_τ_mod[mctr,iz] + dtmp_ϖ[ctr+i] = lin_ϖ_mod[mctr] + dtmp_Z⁺⁺[ctr+i,:,:] .= dAerZ⁺⁺[mctr,:,:] + dtmp_Z⁻⁺[ctr+i,:,:] .= dAerZ⁻⁺[mctr,:,:] + end end - CoreScatteringOpticalProperties(τ_mod, ϖ_mod,AerZ⁺⁺, AerZ⁻⁺), - linCoreScatteringOpticalProperties(lin_τ_mod, lin_ϖ_mod, dAerZ⁺⁺, dAerZ⁻⁺) + + aer = [scatt_paired_xdx(CoreScatteringOpticalProperties( + (τ_mod[i]), ϖ_mod, + AerZ⁺⁺, AerZ⁻⁺), + linCoreScatteringOpticalProperties( + (dtmp_τ[i]), + dtmp_ϖ, + dtmp_Z⁺⁺, dtmp_Z⁻⁺)) + for i=1:Nz] + + #lin_aer = [linCoreScatteringOpticalProperties( + # (dtmp_τ[iz]), + # dtmp_ϖ, + # dtmp_Z⁺⁺, dtmp_Z⁻⁺) for iz=1:Nz] + + return aer #, lin_aer + #CoreScatteringOpticalProperties(τ_mod, ϖ_mod, AerZ⁺⁺, AerZ⁻⁺), + #linCoreScatteringOpticalProperties(dtmp_τ, dtmp_ϖ, dtmp_Z⁺⁺, dtmp_Z⁻⁺) + + #CoreScatteringOpticalProperties(τ_mod, ϖ_mod, AerZ⁺⁺, AerZ⁻⁺), + #linCoreScatteringOpticalProperties(lin_τ_mod, lin_ϖ_mod, dAerZ⁺⁺, dAerZ⁻⁺) end # Extract scattering definitions and integrated absorptions for the source function! @@ -132,24 +275,27 @@ function extractEffectiveProps( #FT = eltype(lods[1].τ) nSpec = length(lods[1].τ) - nZ = length(lods) + Nz = length(lods) + nparams = size(lin_lods[1].lin_τ, 1) # First the Scattering Interfaces: scattering_interface = ScatteringInterface_00() scattering_interfaces_all = [] - τ_sum_all = similar(lods[1].τ,(nSpec,nZ+1)) #?? - #lin_τ_sum_all = similar(lin_lods[1].τ,(nSpec,nZ+1)) #?? + τ_sum_all = similar(lods[1].τ,(nSpec,Nz+1)) #?? + lin_τ_sum_all = similar(lin_lods[1].lin_τ,(nparams,nSpec,Nz+1)) #?? τ_sum_all[:,1] .= 0 - lin_τ_sum_all[:,1] .= 0 + lin_τ_sum_all[:,:,1] .= 0 #@show FT - for iz =1:nZ + for iz =1:Nz # Need to check max entries in Z matrices here as well later! scatter = maximum(lods[iz].τ .* lods[iz].ϖ) > 2eps(FT) - scattering_interface = get_scattering_interface(scattering_interface, scatter, iz) + scattering_interface = get_scattering_interface(scattering_interface, + scatter, iz) push!(scattering_interfaces_all, scattering_interface) #@show typeof(τ_sum_all[:,iz]), typeof(lods[iz].τ) - @views τ_sum_all[:,iz+1] = τ_sum_all[:,iz] + getG_atSun(lods[iz], quad_points) * lods[iz].τ - for ctr = 1:size(lin_τ_sum_all,1) - @views lin_τ_sum_all[ctr,:,iz+1] = lin_τ_sum_all[ctr,:,iz] + getG_atSun(lods[iz], quad_points) * lods[iz].lin_τ[ctr,:] + @views τ_sum_all[:,iz+1] = τ_sum_all[:,iz] + + getG_atSun(lods[iz], quad_points) * lods[iz].τ + for ctr = 1:nparams + @views lin_τ_sum_all[ctr,:,iz+1] = lin_τ_sum_all[ctr,:,iz] #+ getG_atSun(lods[iz], quad_points) * lods[iz].lin_τ[ctr,:] end end return scattering_interfaces_all, τ_sum_all, lin_τ_sum_all @@ -166,7 +312,7 @@ function getG_atSun(lod::CoreDirectionalScatteringOpticalProperties,quad_points: end =# -function (in::linCoreScatteringOpticalProperties, arr_type) +function expandOpticalProperties(in::linCoreScatteringOpticalProperties, arr_type) @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺ = in @assert size(lin_τ) == size(lin_ϖ) "τ and ϖ sizes need to match" if size(lin_Z⁺⁺,4) == 1 @@ -179,18 +325,18 @@ function (in::linCoreScatteringOpticalProperties, arr_type) end end -function expandOpticalProperties(in::linCoreDirectionalScatteringOpticalProperties, arr_type) +function expandOpticalProperties_lin(in::linCoreDirectionalScatteringOpticalProperties, arr_type) @unpack lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺, lin_G = in @assert size(lin_τ) == size(lin_ϖ) "τ and ϖ sizes need to match" - if size(lin_Z⁺⁺,4) == 1 - lin_Z⁺⁺ = _repeat(lin_Z⁺⁺,1,1,1,length(lin_τ[1,:])) - lin_Z⁻⁺ = _repeat(lin_Z⁻⁺,1,1,1,length(lin_τ[1,:])) + if size(lin_Z⁺⁺,5) == 1 + lin_Z⁺⁺ = _repeat(lin_Z⁺⁺,1,1,1,1,length(lin_τ[1,:])) + lin_Z⁻⁺ = _repeat(lin_Z⁻⁺,1,1,1,1,length(lin_τ[1,:])) return linCoreDirectionalScatteringOpticalProperties( arr_type(lin_τ), arr_type(lin_ϖ), arr_type(lin_Z⁺⁺), arr_type(lin_Z⁻⁺), arr_type(lin_G)) else - @assert size(lin_Z⁺⁺,4) == length(lin_τ[1,:]) "Z and τ dimensions need to match " + @assert size(lin_Z⁺⁺,5) == length(lin_τ[1,:]) "Z and τ dimensions need to match " linCoreDirectionalScatteringOpticalProperties( arr_type(lin_τ), arr_type(lin_ϖ), arr_type(lin_Z⁺⁺), arr_type(lin_Z⁻⁺), diff --git a/src/CoreRT/atmo_prof.jl b/src/CoreRT/atmo_prof.jl index 6b5f9bbc..4de598cf 100644 --- a/src/CoreRT/atmo_prof.jl +++ b/src/CoreRT/atmo_prof.jl @@ -5,7 +5,9 @@ This file contains functions that are related to atmospheric profile calculation =# "Compute pressure levels, vmr, vcd for atmospheric profile, given p_half, T, q" -function compute_atmos_profile_fields(T::AbstractArray{FT,1}, p_half::AbstractArray{FT,1}, q, vmr; g₀=9.807) where FT +function compute_atmos_profile_fields(T::AbstractArray{FT,1}, + p_half::AbstractArray{FT,1}, + q, vmr; g₀=9.807) where FT #@show "Atmos", FT # Floating type to use #FT = eltype(T) @@ -15,7 +17,7 @@ function compute_atmos_profile_fields(T::AbstractArray{FT,1}, p_half::AbstractAr p_full = (p_half[2:end] + p_half[1:end-1]) / 2 # Dry and wet mass - dry_mass = FT(28.9644e-3) # in kg/molec, weighted average for N2 and O2 + dry_mass = FT(28.9644e-3) # in kg/mol, weighted average for N2 and O2 wet_mass = FT(18.01534e-3) # just H2O n_layers = length(T) @@ -188,6 +190,7 @@ function getRayleighLayerOptProp(psurf::FT, λ::Union{Array{FT}, FT}, depol_fct: return τRayl end + """ $(FUNCTIONNAME)(total_τ, p₀, σp, p_half) @@ -224,6 +227,7 @@ function getAerosolLayerOptProp(total_τ::FT, dist::Distribution, profile::Atmos τAer = (total_τ / sum(ρ)) * ρ end + """ $(FUNCTIONNAME)(τRayl, τAer, aerosol_optics, Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, τ_abs, arr_type) @@ -304,6 +308,27 @@ function construct_atm_layer(τRayl, τAer, τ_λ = τ_abs .+ τ ϖ_λ = (τ * ϖ) ./ τ_λ + for i = 1:nAer + #@show τ, ϖ , A, τAer[i] + τ += τAer[i] + ϖ += τAer[i] * aerosol_optics[i].ω̃ + A += τAer[i] * aerosol_optics[i].ω̃ * (1 - aerosol_optics[i].fᵗ) + Z⁺⁺ += τAer[i] * aerosol_optics[i].ω̃ * (1 - aerosol_optics[i].fᵗ) * Aer𝐙⁺⁺[:,:,i] + Z⁻⁺ += τAer[i] * aerosol_optics[i].ω̃ * (1 - aerosol_optics[i].fᵗ) * Aer𝐙⁻⁺[:,:,i] + #@show τ, ϖ , A + #for j=1:4 + # ctr = j + (nAer-1)*4 + # if(i==1) + # tmp_lin_τ[ctr] = 1 + # tmp_lin_ϖ[ctr] = τAer[i] + # tmp_lin_A[ctr] = + #end + #for j = 1:7 + # ctr = j + (nAer-1)*7 + # xlin_τ[ctr] = ... + + end + return Array(τ_λ), Array(ϖ_λ), τ, ϖ, Array(Z⁺⁺), Array(Z⁻⁺), fscattRayl end @@ -327,13 +352,13 @@ function construct_all_atm_layers( Z⁺⁺_all = zeros(FT_phase, NquadN, NquadN, Nz) Z⁻⁺_all = zeros(FT_phase, NquadN, NquadN, Nz) - dτ_max_all = zeros(FT_ext, Nz) - dτ_all = zeros(FT_ext, Nz) - fscattRayl_all = zeros(FT_ext, Nz) - ndoubl_all = zeros(Int64, Nz) - dτ_λ_all = zeros(FT_ext, nSpec, Nz) - expk_all = zeros(FT_ext, nSpec, Nz) - scatter_all = zeros(Bool, Nz) + #dτ_max_all = zeros(FT_ext, Nz) + #dτ_all = zeros(FT_ext, Nz) + #fscattRayl_all = zeros(FT_ext, Nz) + #ndoubl_all = zeros(Int64, Nz) + #dτ_λ_all = zeros(FT_ext, nSpec, Nz) + #expk_all = zeros(FT_ext, nSpec, Nz) + #scatter_all = zeros(Bool, Nz) for iz=1:Nz diff --git a/src/CoreRT/lin_atmo_prof.jl b/src/CoreRT/lin_atmo_prof.jl new file mode 100644 index 00000000..4aade944 --- /dev/null +++ b/src/CoreRT/lin_atmo_prof.jl @@ -0,0 +1,486 @@ +#= + +This file contains functions that are related to atmospheric profile calculations + +=# + +"Compute pressure levels, vmr, vcd for atmospheric profile, given p_half, T, q" +function lin_compute_atmos_profile_fields( + T::AbstractArray{FT,1}, + p_half::AbstractArray{FT,1}, + q, vmr,#, + x;#, + #dVMR_CO2, + #dVMR_H2O; + g₀=9.807) where FT + #@show "Atmos", FT + # Floating type to use + #FT = eltype(T) + Nₐ = FT(6.02214179e+23) + R = FT(8.3144598) + + # Calculate full pressure levels + p_full = (p_half[2:end] + p_half[1:end-1]) / 2 + + # Dry and wet mass + dry_mass = FT(28.9644e-3) # in kg/molec, weighted average for N2 and O2 + wet_mass = FT(18.01534e-3) # just H2O + n_layers = length(T) + + # Also get a VMR vector of H2O (volumetric!) + vmr_h2o = zeros(FT, n_layers, ) + vcd_dry = zeros(FT, n_layers, ) + vcd_h2o = zeros(FT, n_layers, ) + Δz = zeros(FT, n_layers) + z = zeros(FT, n_layers) + + psurf = x[1] + @assert x[1]==p_half[end] + # Now actually compute the layer VCDs + M = FT(0.0) + for i = 1:n_layers + Δp = p_half[i + 1] - p_half[i] + a = (i<=65) ? x[8] : x[9] + vmr_h2o[i] = a*dry_mass/(dry_mass-wet_mass*(1-1/q[i])) + vmr_dry = 1 - vmr_h2o[i] + M = vmr_dry * dry_mass + vmr_h2o[i] * wet_mass + vcd = Nₐ * Δp / (M * g₀ * 100^2) * 100 + vcd_dry[i] = vmr_dry * vcd # includes m2->cm2 + vcd_h2o[i] = vmr_h2o[i] * vcd + Δz[i] = (log(p_half[i + 1]/p_half[i])) / (g₀ * M / (R * T[i]) ) + z[1:i] = z[1:i] .+ Δz[i]#@show Δz, T[i], M, Δp + end + #Δp_surf = p_half[end] - p_half[end-1] + dΔz0dpsurf = (1 ./ (p_half[end])) / (g₀ * M / (R * T[end]) ) + dzdpsurf = zeros(length(z)) .+ dΔz0dpsurf + + prof = LogNormal(x[6], x[7]) + vmr["CO2"] = (x[2].+zeros(length(z))) + + (x[3] * exp.(-z./x[5])) + + (x[4] * pdf.(prof, z)) + vmr_co2 = vmr["CO2"] + #dVMR_H2O[1,:] = 0.0 + #dVMR_H2O[1,end] = dVMR_H2O[end]./Δp_surf + dVMR_H2O = zeros(2, length(z)) + dVMR_CO2 = zeros(7, length(z)) + dVMR_H2O[1,:] = [vcd_h2o[1:65]/x[8]; vcd_h2o[66:end] * 0.0;] # wrt x[7] + dVMR_H2O[2,:] = [vcd_h2o[1:65] * 0.0; vcd_h2o[66:end]/x[9];] # wrt x[8] + + dVMR_CO2[1,:] = (x[3] * exp.(-z./x[5]) * (-1/x[5]) .- + (pdf.(prof,z)./z) .* (1 .+ log.(z)/x[7]^2)) .* dzdpsurf; # wrt x[1] + dVMR_CO2[2,:] = 1.0 .+ zeros(length(z)) # wrt x[2] + dVMR_CO2[3,:] = exp.(-z./x[5]) # wrt x[3] + dVMR_CO2[4,:] = pdf.(prof, z) # wrt x[4] + dVMR_CO2[5,:] = x[3] * exp.(-z./x[5]) .* z./(x[5])^2 # wrt x[5] + dVMR_CO2[6,:] = x[4] * pdf.(prof, z) .* (log.(z) .- x[6]) / x[7]^2 + dVMR_CO2[7,:] = (x[4] * pdf.(prof, z) / x[7]) .* + (((log.(z) .- x[6]) / x[7]).^2 .- 1) + + + #= + # TODO: This is still a bit clumsy: + new_vmr = Dict{String, Union{Real, Vector}}() + + for molec_i in keys(vmr) + if vmr[molec_i] isa AbstractArray + + pressure_grid = collect(range(minimum(p_full), maximum(p_full), length=length(vmr[molec_i]))) + interp_linear = LinearInterpolation(pressure_grid, vmr[molec_i]) + new_vmr[molec_i] = [interp_linear(x) for x in p_full] + else + new_vmr[molec_i] = vmr[molec_i] + end + end + =# + #return p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr, Δz, z + return p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, vmr_co2, Δz, z, dzdpsurf, dVMR_H2O, dVMR_CO2 +end +#= +"From a yaml file, get the stored fields (psurf, T, q, ak, bk), calculate derived fields, +and return an AtmosphericProfile object" +function read_atmos_profile(file_path::String) + + # Make sure file is yaml type + @assert endswith(file_path, ".yaml") "File must be yaml" + + # Read in the data and pass to compute fields + params_dict = YAML.load_file(file_path) + + # Validate the parameters before doing anything else + # validate_atmos_profile(params_dict) + + T = convert.(Float64, params_dict["T"]) + + # Calculate derived fields + if ("ak" in keys(params_dict)) + psurf = convert(Float64, params_dict["p_surf"]) + q = convert.(Float64, params_dict["q"]) + ak = convert.(Float64, params_dict["ak"]) + bk = convert.(Float64, params_dict["bk"]) + p_half = (ak + bk * psurf) + p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, Δz = compute_atmos_profile_fields(T, p_half, q, Dict()) + elseif ("q" in keys(params_dict)) + p_half = convert(Float64, params_dict["p_half"]) + psurf = p_half[end] + q = convert.(Float64, params_dict["q"]) + p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, Δz = compute_atmos_profile_fields(T, p_half, q, Dict()) + else + p_half = convert.(Float64, params_dict["p_half"]) + psurf = p_half[end] + q = zeros(length(T)) + p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, Δz = compute_atmos_profile_fields(T, p_half, q, Dict()) + end + + # Convert vmr to appropriate type + vmr = convert(Dict{String, Union{Real, Vector}}, params_dict["vmr"]) + + # Return the atmospheric profile struct + return AtmosphericProfile(T, q, p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, vmr) + +end +=# +"Reduce profile dimensions by re-averaging to near-equidistant pressure grid" +function lin_reduce_profile(n::Int, linprofile::linAtmosphericProfile{FT}) where {FT} + + # Can only reduce the profile, not expand it + @assert n < length(linprofile.T) + + # Unpack the profile vmr + #@unpack vmr, Δz = linprofile + @unpack Δz = linprofile + + # New rough half levels (boundary points) + a = range(0, maximum(linprofile.p_half), length=n+1) + + # Matrices to hold new values + T = zeros(FT, n); + q = zeros(FT, n); + Δz_ = zeros(FT, n); + p_full = zeros(FT, n); + p_half = zeros(FT, n+1); + z = zeros(FT, n); + vmr_h2o = zeros(FT, n); + vmr_co2 = zeros(FT, n); + vcd_dry = zeros(FT, n); + vcd_h2o = zeros(FT, n); + dzdpsurf = zeros(FT, n); + dVMR_H2O = zeros(FT, 2, n); + dVMR_CO2 = zeros(FT, 7, n); + + # Loop over target number of layers + indices = [] + for i = 1:n + + # Get the section of the atmosphere with the i'th section pressure values + ind = findall(a[i] .< linprofile.p_full .<= a[i+1]); + push!(indices, ind) + @assert length(ind) > 0 "Profile reduction has an empty layer" + #@show i, ind, a[i], a[i+1] + # Set the pressure levels accordingly + p_half[i] = a[i] # profile.p_half[ind[1]] + p_half[i+1] = a[i+1] # profile.p_half[ind[end]] + + # Re-average the other parameters to produce new layers + + + p_full[i] = mean(linprofile.p_full[ind]) + T[i] = mean(linprofile.T[ind]) + q[i] = mean(linprofile.q[ind]) + Δz_[i] = sum(Δz[ind]) + z[i] = maximum(linprofile.z[ind]) + vcd_dry[i] = sum(linprofile.vcd_dry[ind]) + vcd_h2o[i] = sum(linprofile.vcd_h2o[ind]) + vmr_h2o[i] = sum(linprofile.vmr_h2o[ind].*linprofile.p_half[ind]./linprofile.T[ind])/ + sum(linprofile.p_half[ind]./linprofile.T[ind])#vcd_h2o[i]/vcd_dry[i] + vmr_co2[i] = sum(linprofile.vmr_co2[ind].*linprofile.p_half[ind]./linprofile.T[ind])/ + sum(linprofile.p_half[ind]./linprofile.T[ind]) + dzdpsurf[i] = mean(linprofile.dzdpsurf[ind]) + for j=1:2 + dVMR_H2O[j,i] = sum(linprofile.dVMR_H2O[j,ind].*linprofile.p_half[ind]./linprofile.T[ind])/ + sum(linprofile.p_half[ind]./linprofile.T[ind]) + dVMR_CO2[j,i] = sum(linprofile.dVMR_CO2[j,ind].*linprofile.p_half[ind]./linprofile.T[ind])/ + sum(linprofile.p_half[ind]./linprofile.T[ind]) + end + for j=3:7 + dVMR_CO2[j,i] = sum(linprofile.dVMR_CO2[j,ind].*linprofile.p_half[ind]./linprofile.T[ind])/ + sum(linprofile.p_half[ind]./linprofile.T[ind]) + end + end + #@show indices +#= + new_vmr = Dict{String, Union{Real, Vector}}() + + # need to double check this logic, maybe better to add VCDs?! + for molec_i in keys(vmr) + if profile.vmr[molec_i] isa AbstractArray + # TODO: This needs a VCD_dry weighted average! + new_vmr[molec_i] = [mean(profile.vmr[molec_i][ind]) for ind in indices] + else + new_vmr[molec_i] = profile.vmr[molec_i] + end + end +=# + return linAtmosphericProfile(T, p_full, q, p_half, vmr_h2o, vcd_dry, vcd_h2o, vmr_co2, Δz_, z, dzdpsurf, dVMR_H2O, dVMR_CO2) +end + +""" + $(FUNCTIONNAME)(psurf, λ, depol_fct, vcd_dry) + +Returns the Rayleigh optical thickness per layer at reference wavelength `λ` (N₂,O₂ atmosphere, i.e. terrestrial) + +Input: + - `psurf` surface pressure in `[hPa]` + - `λ` wavelength in `[μm]` + - `depol_fct` depolarization factor + - `vcd_dry` dry vertical column (no water) per layer +""" +function getRayleighLayerOptProp_lin(psurf::FT, λ::Union{Array{FT}, FT}, depol_fct::FT, vcd_dry::Array{FT}) where FT + # TODO: Use noRS/noRS_plus to use n2/o2 molecular constants + # to compute tau_scat and depol_fct + Nz = length(vcd_dry) + τRayl = zeros(FT,size(λ,1),Nz) + lin_τRayl = zeros(FT,size(λ,1),Nz) # derivative of τRayl wrt psurf + # Total vertical Rayleigh scattering optical thickness, TODO: enable sub-layers and use VCD based taus + tau_scat = FT(0.00864) * (psurf / FT(1013.25)) * λ.^(-FT(3.916) .- FT(0.074) * λ .- FT(0.05) ./ λ) + tau_scat = tau_scat * (FT(6.0) + FT(3.0) * depol_fct) / (FT(6.0)- FT(7.0) * depol_fct) + # @show tau_scat, λ + k = tau_scat / sum(vcd_dry) + for i = 1:Nz + τRayl[:,i] .= k * vcd_dry[i] + lin_τRayl[:,i] .= τRayl[:,i]/psurf + end + return τRayl, lin_τRayl +end + + +""" + $(FUNCTIONNAME)(total_τ, p₀, σp, p_half) + +Returns the aerosol optical depths per layer using a Gaussian distribution function with p₀ and σp on a pressure grid +""" + +function getAerosolLayerOptProp_lin(total_τ, z₀, σz, z, dzdpsurf)#, p_half) + + # Need to make sure we can also differentiate wrt σp (FT can be Dual!) + FT = eltype(z₀) + Nz = length(z) + #ρ = zeros(FT,Nz) + #dρdz₀ = zeros(FT,Nz) + #dρdσz = zeros(FT,Nz) + # @show p_half, p₀, σp + + prof = LogNormal(log(z₀), σz) + τ_aer = total_τ * pdf.(prof, z) + lin_τ_aer_psurf = - τ_aer./z .* + (1 .+ log.(z)/σz^2) .* dzdpsurf + lin_τ_aer_z₀ = τ_aer .* (log.(z) .- log(z₀)) / σz^2 + lin_τ_aer_σz = (τ_aer / σz) .* + (((log.(z) .- log(z₀)) / σz).^2 .- 1) + + # return convert(FT, τ_aer, lin_τ_aer_psurf, lin_τ_aer_z₀, lin_τ_aer_σz) + return τ_aer, lin_τ_aer_psurf, lin_τ_aer_z₀, lin_τ_aer_σz; + +end + + + +#= +""" + $(FUNCTIONNAME)(τRayl, τAer, aerosol_optics, Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, τ_abs, arr_type) + +Computes the composite layer single scattering parameters (τ, ϖ, Z⁺⁺, Z⁻⁺) + +Returns: + - `τ`, `ϖ` : only Rayleigh scattering and aerosol extinction, no gaseous absorption (no wavelength dependence) + - `τ_λ`,`ϖ_λ`: Rayleigh scattering + aerosol extinction + gaseous absorption (wavelength dependent) + - `Z⁺⁺`,`Z⁻⁺`: Composite Phase matrix (weighted average of Rayleigh and aerosols) + - `fscattRayl`: Rayleigh scattering fraction (needed for Raman computations) +Arguments: + - `τRay` layer optical depth for Rayleigh + - `τAer` layer optical depth for Aerosol(s) (vector) + - `aerosol_optics` array of aerosol optics struct + - `Rayl𝐙⁺⁺` Rayleigh 𝐙⁺⁺ phase matrix (2D) + - `Rayl𝐙⁻⁺` Rayleigh 𝐙⁻⁺ phase matrix (2D) + - `Aer𝐙⁺⁺` Aerosol 𝐙⁺⁺ phase matrix (3D) + - `Aer𝐙⁻⁺` Aerosol 𝐙⁻⁺ phase matrix (3D) + - `τ_abs` layer absorption optical depth array (per wavelength) by gaseous absorption +""" +function construct_atm_layer(τRayl, τAer, + ϖ_Cabannes, #elastic fraction of Rayleigh scattering + aerosol_optics, + Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, + Aer𝐙⁺⁺, Aer𝐙⁻⁺, + τ_abs, arr_type) + + FT = eltype(τRayl) + nAer = length(aerosol_optics) + + # Fixes Rayleigh SSA to 1 for purely elastic (RS_type = noRS) scattering, + # and assumes values less than 1 for Raman scattering + ϖRayl = ϖ_Cabannes #FT(1) + @show ϖRayl + @assert length(τAer) == nAer "Sizes don't match" + + τ = FT(0) + ϖ = FT(0) + A = FT(0) + Z⁺⁺ = similar(Rayl𝐙⁺⁺); + Z⁻⁺ = similar(Rayl𝐙⁺⁺); + + if (τRayl + sum(τAer)) < eps(FT) + fill!(Z⁺⁺, 0); fill!(Z⁻⁺, 0); + return FT(0), FT(1), Z⁺⁺, Z⁻⁺ + end + + τ += τRayl + #@show τRayl, ϖRayl[1], ϖ + ϖ += τRayl * ϖRayl[1] + A += τRayl * ϖRayl[1] + + Z⁺⁺ = τRayl * ϖRayl[1] * Rayl𝐙⁺⁺ + Z⁻⁺ = τRayl * ϖRayl[1] * Rayl𝐙⁻⁺ + + for i = 1:nAer + #@show τ, ϖ , A, τAer[i] + τ += τAer[i] + ϖ += τAer[i] * aerosol_optics[i].ω̃ + A += τAer[i] * aerosol_optics[i].ω̃ * (1 - aerosol_optics[i].fᵗ) + Z⁺⁺ += τAer[i] * aerosol_optics[i].ω̃ * (1 - aerosol_optics[i].fᵗ) * Aer𝐙⁺⁺[:,:,i] + Z⁻⁺ += τAer[i] * aerosol_optics[i].ω̃ * (1 - aerosol_optics[i].fᵗ) * Aer𝐙⁻⁺[:,:,i] + #@show τ, ϖ , A + end + + Z⁺⁺ /= A + Z⁻⁺ /= A + A /= ϖ + ϖ /= τ + + # Rescaling composite SSPs according to Eqs. A.3 of Sanghavi et al. (2013) or Eqs.(8) of Sanghavi & Stephens (2015) + #@show τRayl, τ,A, ϖ + τ *= (FT(1) - (FT(1) - A) * ϖ) + ϖ *= A / (FT(1) - (FT(1) - A) * ϖ)#Suniti + #@show τRayl, τ + fscattRayl = τRayl/τ + # Adding absorption optical depth / albedo: + τ_λ = τ_abs .+ τ + ϖ_λ = (τ * ϖ) ./ τ_λ + + return Array(τ_λ), Array(ϖ_λ), τ, ϖ, Array(Z⁺⁺), Array(Z⁻⁺), fscattRayl +end + +"When performing RT_run, this function pre-calculates properties for all layers, before any Core RT is performed" +function construct_all_atm_layers( + FT, nSpec, Nz, NquadN, + τRayl, τAer, aerosol_optics, + Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, + τ_abs, + ϖ_Cabannes, + arr_type, qp_μ, μ₀, m) + + FT_ext = eltype(τAer) + FT_phase = eltype(τAer) + + # Empty matrices to hold all values + τ_λ_all = zeros(FT_ext, nSpec, Nz) + ϖ_λ_all = zeros(FT_ext, nSpec, Nz) + τ_all = zeros(FT_ext, Nz) + ϖ_all = zeros(FT_ext, Nz) + Z⁺⁺_all = zeros(FT_phase, NquadN, NquadN, Nz) + Z⁻⁺_all = zeros(FT_phase, NquadN, NquadN, Nz) + + dτ_max_all = zeros(FT_ext, Nz) + dτ_all = zeros(FT_ext, Nz) + fscattRayl_all = zeros(FT_ext, Nz) + ndoubl_all = zeros(Int64, Nz) + dτ_λ_all = zeros(FT_ext, nSpec, Nz) + expk_all = zeros(FT_ext, nSpec, Nz) + scatter_all = zeros(Bool, Nz) + + for iz=1:Nz + + # Construct atmospheric properties + τ_λ_all[:, iz], + ϖ_λ_all[:, iz], + τ_all[iz], + ϖ_all[iz], + Z⁺⁺_all[:,:,iz], + Z⁻⁺_all[:,:,iz], + fscattRayl_all[iz] = construct_atm_layer(τRayl[iz], τAer[:,iz], + ϖ_Cabannes, + aerosol_optics, + Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, + τ_abs[:,iz], arr_type) + #@show fscattRayl_all[iz] + # Compute doubling number + dτ_max_all[iz] = minimum([τ_all[iz] * ϖ_all[iz], FT(0.001) * minimum(qp_μ)]) + dτ_all[iz], ndoubl_all[iz] = doubling_number(dτ_max_all[iz], τ_all[iz] * ϖ_all[iz]) #Suniti + + # Compute dτ vector + dτ_λ_all[:, iz] = (τ_λ_all[:, iz] ./ (FT(2)^ndoubl_all[iz])) + #@show maximum(dτ_λ_all[:,iz]) + expk_all[:, iz] = exp.(-dτ_λ_all[:, iz] /μ₀) #Suniti + + # Determine whether there is scattering + scatter_all[iz] = ( sum(τAer[:,iz]) > 1.e-8 || + (( τRayl[iz] > 1.e-8 ) && (m < 3))) ? + true : false + end + + # Compute sum of optical thicknesses of all layers above the current layer + τ_sum_all = accumulate(+, τ_λ_all, dims=2) + + # First start with all zeros + # At the bottom of the atmosphere, we have to compute total τ_sum (bottom of lowest layer), for the surface interaction + τ_sum_all = hcat(zeros(FT, size(τ_sum_all[:,1])), τ_sum_all) + + # Starting scattering interface (None for both added and composite) + scattering_interface = ScatteringInterface_00() + scattering_interfaces_all = [] + + for iz = 1:Nz + # Whether there is scattering in the added layer, composite layer, neither or both + scattering_interface = get_scattering_interface(scattering_interface, scatter_all[iz], iz) + push!(scattering_interfaces_all, scattering_interface) + end + + return ComputedAtmosphereProperties(τ_λ_all, ϖ_λ_all, τ_all, ϖ_all, Z⁺⁺_all, Z⁻⁺_all, dτ_max_all, dτ_all, ndoubl_all, dτ_λ_all, expk_all, scatter_all, τ_sum_all, fscattRayl_all, scattering_interfaces_all) +end +=# + +# TODO: +"Given the CrossSectionModel, the grid, and the AtmosphericProfile, fill up the τ_abs array with the cross section at each layer +(using pressures/temperatures) from the profile" +function compute_absorption_profile_lin!(τ_abs::Array{FT,2}, + lin_τ_abs::Array{FT,3}, + Δp_surf, + dVMR, + #dVMR_CO2, + absorption_model, + grid, + vmr, + profile::linAtmosphericProfile, + ) where FT + + # The array to store the cross-sections must be same length as number of layers + @assert size(τ_abs,2) == length(profile.p_full) + @assert length(vmr) ==1 || length(vmr) == length(profile.p_full) "Length of VMR array has to match profile size or be uniform" + #@show grid + @showprogress 1 for iz in 1:length(profile.p_full) + + # Pa -> hPa + p = profile.p_full[iz] + T = profile.T[iz] + # Either use the current layer's vmr, or use the uniform vmr + vmr_curr = vmr isa AbstractArray ? vmr[iz] : vmr + Δτ = Array(absorption_cross_section(absorption_model, grid, p, T)) * profile.vcd_dry[iz] * vmr_curr + τ_abs[:,iz] += Δτ # Array(absorption_cross_section(absorption_model, grid, p, T)) * profile.vcd_dry[iz] * vmr_curr + + for ipar in 1:9 + lin_τ_abs[ipar,:,iz] += Δτ * (dVMR[ipar,iz]./vmr_curr) + end + if iz==length(profile.p_full) + lin_τ_abs[1,:,iz] += Δτ/Δp_surf + end + end + +end \ No newline at end of file diff --git a/src/CoreRT/lin_model_from_parameters.jl b/src/CoreRT/lin_model_from_parameters.jl index 81f155ab..f5eea3d6 100644 --- a/src/CoreRT/lin_model_from_parameters.jl +++ b/src/CoreRT/lin_model_from_parameters.jl @@ -9,7 +9,10 @@ like optical thicknesses, from the input parameters. Produces a vSmartMOM_Model default_parameters() = parameters_from_yaml(joinpath(dirname(pathof(vSmartMOM)), "CoreRT", "DefaultParameters.yaml")) "Take the parameters specified in the vSmartMOM_Parameters struct, and calculate derived attributes into a vSmartMOM_Model" -function lin_model_from_parameters(params::vSmartMOM_Parameters) +function lin_model_from_parameters(params::vSmartMOM_Parameters, + x)#, + #dVMR_CO2, + #dVMR_H2O) FT = params.float_type #@show FT # Number of total bands and aerosols (for convenience) @@ -27,26 +30,47 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) # Get AtmosphericProfile from parameters vmr = isnothing(params.absorption_params) ? Dict() : params.absorption_params.vmr - p_full, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr, Δz = compute_atmos_profile_fields(params.T, params.p, params.q, vmr) - - profile = AtmosphericProfile(params.T, p_full, params.q, p_half, vmr_h2o, vcd_dry, vcd_h2o, new_vmr,Δz) + p_full, p_half, vmr_h2o, vcd_dry, + vcd_h2o, vmr_co2, Δz, z, dzdpsurf, dVMR_H2O, dVMR_CO2 = + lin_compute_atmos_profile_fields(params.T, + params.p, + params.q, + vmr, + x)#, + #dVMR_CO2, + #dVMR_H2O) + + + linprofile = linAtmosphericProfile(params.T, + p_full, + params.q, + p_half, + vmr_h2o, + vcd_dry, vcd_h2o, + vmr_co2, + Δz, z, dzdpsurf, + dVMR_H2O, dVMR_CO2) # Reduce the profile to the number of target layers (if specified) if params.profile_reduction_n != -1 - profile = reduce_profile(params.profile_reduction_n, profile); + linprofile = lin_reduce_profile(params.profile_reduction_n, linprofile); end # Rayleigh optical properties calculation greek_rayleigh = Scattering.get_greek_rayleigh(FT(params.depol)) - # Remove rayleight for testing: - τ_rayl = [zeros(FT,length(params.spec_bands[i]), length(params.T)) for i=1:n_bands]; + # Remove rayleigh for testing: + τ_rayl = [zeros(FT,length(params.spec_bands[i]), length(linprofile.T)) for i=1:n_bands]; + lin_τ_rayl = [zeros(FT,length(params.spec_bands[i]), length(linprofile.T)) for i=1:n_bands]; + #τ_rayl = [zeros(FT,1,length(profile.T)) for i=1:n_bands]; # This is a kludge for now, tau_abs sometimes needs to be a dual. Suniti & us need to rethink this all!! # i.e. code the rt core with fixed amount of derivatives as in her paper, then compute chain rule for dtau/dVMr, etc... - FT2 = isnothing(params.absorption_params) || !haskey(params.absorption_params.vmr,"CO2") ? params.float_type : eltype(params.absorption_params.vmr["CO2"]) - τ_abs = [zeros(FT2, length(params.spec_bands[i]), length(profile.p_full)) for i in 1:n_bands] - + #FT2 = isnothing(params.absorption_params) || !haskey(params.absorption_params.vmr,"CO2") ? params.float_type : eltype(params.absorption_params.vmr["CO2"]) + τ_abs = [zeros(FT, length(params.spec_bands[i]), length(linprofile.p_full)) for i in 1:n_bands] + lin_τ_abs = [zeros(FT, 9, length(params.spec_bands[i]), length(linprofile.p_full)) for i in 1:n_bands] + VMR = zeros(FT, length(linprofile.T)); + dVMR = zeros(FT, 9, length(linprofile.T)); # Loop over all bands: for i_band=1:n_bands @@ -54,14 +78,15 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) curr_band_λ = FT.(1e4 ./ params.spec_bands[i_band]) # @show profile.vcd_dry, size(τ_rayl[i_band]) # Compute Rayleigh properties per layer for `i_band` band center - - τ_rayl[i_band] .= getRayleighLayerOptProp(profile.p_half[end], + + τ_rayl[i_band], lin_τ_rayl[i_band] = getRayleighLayerOptProp_lin( + linprofile.p_half[end], curr_band_λ, #(mean(curr_band_λ)), - params.depol, profile.vcd_dry); + params.depol, linprofile.vcd_dry); #@show τ_rayl[i_band] # If no absorption, continue to next band isnothing(params.absorption_params) && continue - + Δp_surf = linprofile.p_half[end] - linprofile.p_half[end-1]; # Loop over all molecules in this band, obtain profile for each, and add them up for molec_i in 1:length(params.absorption_params.molecules[i_band]) @show params.absorption_params.molecules[i_band][molec_i] @@ -79,11 +104,40 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) architecture = params.architecture, vmr = 0);#mean(profile.vmr[params.absorption_params.molecules[i_band][molec_i]])) # Calculate absorption profile - - @timeit "Absorption Coeff" compute_absorption_profile!(τ_abs[i_band], absorption_model, params.spec_bands[i_band],profile.vmr[params.absorption_params.molecules[i_band][molec_i]], profile); + if params.absorption_params.modecules[i_band]=="O2" + VMR = linprofile.vmr[params.absorption_params.molecules[i_band][molec_i]] + dVMR .= 0 + elseif params.absorption_params.modecules[i_band]=="CO2" + VMR = linprofile.vmr_co2 + dVMR[1:7,:] .= linprofile.dVMR_CO2 + dVMR[8:9,:] .= 0 + elseif params.absorption_params.modecules[i_band]=="H2O" + VMR = linprofile.vmr_h2o + dVMR[1:7,:] .= 0 + dVMR[8:9,:] .= linprofile.dVMR_H2O + end + @timeit "Absorption Coeff" compute_absorption_profile_lin!( + τ_abs[i_band], + lin_τ_abs[i_band], + Δp_surf, + dVMR, + absorption_model, + params.spec_bands[i_band], + VMR, + #profile.vmr[params.absorption_params.molecules[i_band][molec_i]], + linprofile); # Use LUT directly else - compute_absorption_profile!(τ_abs[i_band], params.absorption_params.luts[i_band][molec_i], params.spec_bands[i_band],profile.vmr[params.absorption_params.molecules[i_band][molec_i]], profile); + compute_absorption_profile_lin!( + τ_abs[i_band], + lin_τ_abs[i_band], + Δp_surf, + dVMR, + params.absorption_params.luts[i_band][molec_i], + params.spec_bands[i_band], + VMR, + #profile.vmr[params.absorption_params.molecules[i_band][molec_i]], + linprofile); end end end @@ -97,11 +151,13 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) #FT2 = params.float_type # τ_aer[iBand][iAer,iZ] - τ_aer = [zeros(FT2, n_aer, length(profile.p_full)) for i=1:n_bands]; + τ_aer = [zeros(FT2, n_aer, length(linprofile.p_full)) for i=1:n_bands]; + lin_τ_aer_psurf = [zeros(FT2, n_aer, length(linprofile.p_full)) for i=1:n_bands]; + lin_τ_aer_z₀ = [zeros(FT2, n_aer, length(linprofile.p_full)) for i=1:n_bands]; + lin_τ_aer_σz = [zeros(FT2, n_aer, length(linprofile.p_full)) for i=1:n_bands]; # Loop over aerosol type for i_aer=1:n_aer - # Get curr_aerosol c_aero = params.scattering_params.rt_aerosols[i_aer] curr_aerosol = c_aero.aerosol @@ -125,7 +181,8 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) mie_model.aerosol.nᵣ = real(params.scattering_params.n_ref) mie_model.aerosol.nᵢ = -imag(params.scattering_params.n_ref) @show params.scattering_params.n_ref - k_ref, dk_ref = compute_ref_aerosol_extinction(mie_model, params.float_type) + k_ref, dk_ref = compute_ref_aerosol_extinction_lin(mie_model)#, + #params.float_type) @show k_ref #params.scattering_params.rt_aerosols[i_aer].p₀, params.scattering_params.rt_aerosols[i_aer].σp # Loop over bands @@ -148,7 +205,7 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) # Compute raw (not truncated) aerosol optical properties (not needed in RT eventually) #@show FT2, FT @timeit "Mie calc" aerosol_optics_raw, lin_aerosol_optics_raw = - compute_aerosol_optical_properties(mie_model, FT); + compute_aerosol_optical_properties_lin(mie_model, FT); @show aerosol_optics_raw.k aerosol_optics_raw.k_ref = k_ref @show k_ref, dk_ref @@ -156,17 +213,21 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) lin_aerosol_optics_raw.dk_ref = dk_ref # Compute truncated aerosol optical properties (phase function and fᵗ), consistent with Ltrunc: #@show i_aer, i_band - aerosol_optics[i_band][i_aer],lin_aerosol_optics[i_band][i_aer] = Scattering.truncate_phase(truncation_type, + aerosol_optics[i_band][i_aer], + lin_aerosol_optics[i_band][i_aer] = Scattering.truncate_phase(truncation_type, aerosol_optics_raw, lin_aerosol_optics_raw; reportFit=false) #aerosol_optics[i_band][i_aer] = aerosol_optics_raw - + @show lin_aerosol_optics[i_band][i_aer].dω̃ @show aerosol_optics[i_band][i_aer].k #@show aerosol_optics[i_band][i_aer].fᵗ # Compute nAer aerosol optical thickness profiles - τ_aer[i_band][i_aer,:] = - params.scattering_params.rt_aerosols[i_aer].τ_ref * - (aerosol_optics[i_band][i_aer].k/k_ref) * - getAerosolLayerOptProp(1, c_aero.profile, profile) + τ_total = params.scattering_params.rt_aerosols[i_aer].τ_ref * + (aerosol_optics[i_band][i_aer].k/k_ref) + τ_aer[i_band][i_aer,:], + lin_τ_aer_psurf[i_band][i_aer, :], + lin_τ_aer_z₀[i_band][i_aer,:], + lin_τ_aer_σz[i_band][i_aer,:] = + getAerosolLayerOptProp_lin(τ_total, exp(c_aero.profile.μ), c_aero.profile.σ, linprofile.z, linprofile.dzdpsurf) @info "AOD at band $i_band : $(sum(τ_aer[i_band][i_aer,:])), truncation factor = $(aerosol_optics[i_band][i_aer].fᵗ)" end end @@ -196,7 +257,15 @@ function lin_model_from_parameters(params::vSmartMOM_Parameters) τ_rayl, τ_aer, obs_geom, - profile), vSmartMOM_lin(lin_aerosol_optics) + linprofile), + vSmartMOM_lin(lin_aerosol_optics, + lin_τ_rayl, + lin_τ_abs, + lin_τ_aer_psurf, + lin_τ_aer_z₀, + lin_τ_aer_σz) + + # TODO: add dτ_abs[1] due to O2 end diff --git a/src/CoreRT/lin_types.jl b/src/CoreRT/lin_types.jl index 384f39bd..0e2b095b 100644 --- a/src/CoreRT/lin_types.jl +++ b/src/CoreRT/lin_types.jl @@ -15,9 +15,9 @@ This file contains the linearization of all relevant types that are used in the - `ComputedAtmosphereProperties` and `ComputedLayerProperties` hold intermediate computed properties =# -#= + "Struct for an atmospheric profile" -struct AtmosphericProfile{FT, VMR <: Union{Real, Vector}} +struct linAtmosphericProfile{FT}#{FT, VMR <: Union{Real, Vector}} "Temperature Profile" T::Array{FT,1} "Pressure Profile (Full)" @@ -33,11 +33,18 @@ struct AtmosphericProfile{FT, VMR <: Union{Real, Vector}} "Vertical Column Density (H2O)" vcd_h2o::Array{FT,1} "Volume Mixing Ratio of Constituent Gases" - vmr::Dict{String, VMR} + vmr_co2::Array{FT,1} #Dict{String, VMR} "Layer height (meters)" Δz::Array{FT,1} + "Layer altitude (meters)" + z::Array{FT,1} + dzdpsurf::Array{FT,1} + "Derivative of vmr_h2o with respect to multiplicative factors in the boundary layer and above" + dVMR_H2O::Array{FT,2} + "Derivative of vmr_co2 with respect to psurf, 3 multiplicative factors for uniform, exponential, and lognormal profiles, and their respective parameters H, z₀, and σ" + dVMR_CO2::Array{FT,2} end - +#= "Types for describing atmospheric parameters" abstract type AbstractObsGeometry end @@ -510,32 +517,36 @@ struct QuadPoints{FT} end =# """ - struct vSmartMOM_Model + struct vSmartMOM_lin A struct which holds all derived model parameters (including any computations) # Fields $(DocStringExtensions.FIELDS) """ -mutable struct vSmartMOM_lin{linAE}#, linTAB, linTR, linTAE, linPRO} +mutable struct vSmartMOM_lin #<:Union{Vector{Vector{dAerosolOptics}}, Vector{Array{FT},1}, Vector{Array{FT},3}} where FT#{FT} <: AbstractFloat#{linAE, linTAE, linTAE}#, linTAB, linTR, linTAE, linPRO} #"Struct with all individual parameters" #params::PA # vSmartMOM_Parameters "Truncated aerosol optics" - lin_aerosol_optics::linAE # AbstractArray{AbstractArray{AerosolOptics}} + #lin_aerosol_optics::dAerosolOptics{FT} #::linAE # AbstractArray{AbstractArray{AerosolOptics}} + lin_aerosol_optics::Any #Vector{Vector{dAerosolOptics}} #"Greek coefs in Rayleigh calculations" #greek_rayleigh::GR # GreekCoefs #"Quadrature points/weights, etc" #quad_points::QP # QuadPoints #"Array to hold cross-sections over entire atmospheric profile" - #lin_τ_abs::linTAB # AbstractArray{AbstractArray} + #"Rayleigh optical thickness" - #lin_τ_rayl::linTR # AbstractArray{AbstractArray} + lin_τ_rayl::Any #Vector{Array{FT},1} #Array{Matrix{FT},1}#linTR # AbstractArray{AbstractArray} + #"Molecular absorption" + lin_τ_abs::Any #Vector{Array{FT},3} #linTAB # AbstractArray{AbstractArray} #"Aerosol optical thickness" - #lin_τ_aer::linTAE # AbstractArray{AbstractArray} - + lin_τ_aer_psurf::Any #Vector{Array{FT},1} #Array{Matrix{FT},1} + lin_τ_aer_z₀::Any #Vector{Array{FT},1}#Array{Matrix{FT},1}#::linTAE # AbstractArray{AbstractArray} + lin_τ_aer_σz::Any #Vector{Array{FT},1}#Array{Matrix{FT},1}#::linTAE #"Observational Geometry (includes sza, vza, vaz)" #obs_geom::Ogeom # ObsGeometry #"Atmospheric profile to use" @@ -632,10 +643,10 @@ Base.@kwdef struct linComputedLayerProperties #scattering_interface end -abstract type AbstractOpticalProperties end +abstract type linAbstractOpticalProperties end # Core optical Properties COP -Base.@kwdef struct linCoreScatteringOpticalProperties{FT,FT2,FT3} <: AbstractOpticalProperties +Base.@kwdef struct linCoreScatteringOpticalProperties{FT,FT2,FT3} <: linAbstractOpticalProperties "Absorption optical depth (scalar or wavelength dependent)" lin_τ::FT "Single scattering albedo" @@ -647,7 +658,7 @@ Base.@kwdef struct linCoreScatteringOpticalProperties{FT,FT2,FT3} <: AbstractOp end # Core optical Properties COP with directional cross section -Base.@kwdef struct linCoreDirectionalScatteringOpticalProperties{FT,FT2,FT3,FT4} <: AbstractOpticalProperties +Base.@kwdef struct linCoreDirectionalScatteringOpticalProperties{FT,FT2,FT3,FT4} <: linAbstractOpticalProperties "Absorption optical depth (scalar or wavelength dependent)" lin_τ::FT "Single scattering albedo" @@ -660,69 +671,125 @@ Base.@kwdef struct linCoreDirectionalScatteringOpticalProperties{FT,FT2,FT3,FT4} lin_G::FT4 end -Base.@kwdef struct linCoreAbsorptionOpticalProperties{FT} <: AbstractOpticalProperties +Base.@kwdef struct linCoreAbsorptionOpticalProperties{FT} <: linAbstractOpticalProperties "Absorption optical depth (scalar or wavelength dependent)" lin_τ::FT end +struct paired_xdx + x::AbstractOpticalProperties #{xFT, xFT2, xFT3} + dx::linAbstractOpticalProperties #{dxFT, dxFT2, dxFT3} + #x::Vector{AbstractOpticalProperties} #{xFT, xFT2, xFT3} + #dx::Vector{linAbstractOpticalProperties} #{dxFT, dxFT2, dxFT3} +end + +struct scatt_paired_xdx + x::CoreScatteringOpticalProperties #{xFT, xFT2, xFT3} + dx::linCoreScatteringOpticalProperties #{dxFT, dxFT2, dxFT3} + #x::Vector{CoreScatteringOpticalProperties} #{xFT, xFT2, xFT3} + #dx::Vector{linCoreScatteringOpticalProperties} #{dxFT, dxFT2, dxFT3} +end + +struct dirscatt_paired_xdx + x::CoreDirectionalScatteringOpticalProperties #{xFT, xFT2, xFT3} + dx::linCoreDirectionalScatteringOpticalProperties #{dxFT, dxFT2, dxFT3} + #x::Vector{CoreDirectionalScatteringOpticalProperties} #{xFT, xFT2, xFT3} + #dx::Vector{linCoreDirectionalScatteringOpticalProperties} #{dxFT, dxFT2, dxFT3} +end + +struct abs_paired_xdx + x::CoreAbsorptionOpticalProperties #{xFT, xFT2, xFT3} + dx::linCoreAbsorptionOpticalProperties #{dxFT, dxFT2, dxFT3} + #x::Vector{CoreAbsorptionOpticalProperties} #{xFT, xFT2, xFT3} + #dx::Vector{linCoreAbsorptionOpticalProperties} #{dxFT, dxFT2, dxFT3} +end + # Adding Core Optical Properties, can have mixed dimensions! -function Base.:+( x::CoreScatteringOpticalProperties{xFT, xFT2, xFT3}, - y::CoreScatteringOpticalProperties{yFT, yFT2, yFT3}, - dx::linCoreScatteringOpticalProperties{xFT, xFT2, xFT3}, - dy::linCoreScatteringOpticalProperties{yFT, yFT2, yFT3} - ) where {xFT, xFT2, xFT3, yFT, yFT2, yFT3} - # Predefine some arrays: +function Base.:+(z1::scatt_paired_xdx, + z2::scatt_paired_xdx)# where {xFT, xFT2, xFT3, yFT, yFT2, yFT3} + #( (x::CoreScatteringOpticalProperties{xFT, xFT2, xFT3}, + # dx::linCoreScatteringOpticalProperties{dxFT, dxFT2, dxFT3}), + # (y::CoreScatteringOpticalProperties{yFT, yFT2, yFT3}, + # dy::linCoreScatteringOpticalProperties{dyFT, dyFT2, dyFT3}) + # ) where {xFT, xFT2, xFT3, yFT, yFT2, yFT3} + # Predefine some arrays: + x = z1.x + dx = z1.dx + y = z2.x + dy = z2.dx + xZ⁺⁺ = x.Z⁺⁺ xZ⁻⁺ = x.Z⁻⁺ yZ⁺⁺ = y.Z⁺⁺ yZ⁻⁺ = y.Z⁻⁺ τ = x.τ .+ y.τ - wx = x.τ .* x.ϖ + wx = x.τ .* x.ϖ wy = y.τ .* y.ϖ - w = wx .+ wy - ϖ = w ./ τ - - #xlin_Z⁺⁺ = dx.lin_Z⁺⁺ - #xlin_Z⁻⁺ = dx.lin_Z⁻⁺ - #ylin_Z⁺⁺ = dy.lin_Z⁺⁺ - #ylin_Z⁻⁺ = dy.lin_Z⁻⁺ - - lin_τ = [dx.lin_τ, dy.lin_τ] - lin_wx = (dx.lin_τ .* (x.ϖ - ϖ) .+ x.τ .* dx.lin_ϖ)./τ - lin_wy = (dy.lin_τ .* (y.ϖ - ϖ) .+ y.τ .* dy.lin_ϖ)./τ - lin_ϖ = [lin_wx, lin_wy] + @show size(wx), size(wy), size(τ) + ϖ = (wx .+ wy) ./ τ #(wx + wy) ./ τ + + @show size(dx.lin_τ), size(dy.lin_τ) + nparams = size(dx.lin_τ,2) + lin_τ = similar(dx.lin_τ) + lin_ϖ = similar(dx.lin_τ) + #for i = 1:nparams + lin_τ = dx.lin_τ .+ dy.lin_τ #dx.lin_τ + dy.lin_τ + @show size(lin_ϖ), size(dx.lin_τ), size(x.ϖ), size(dx.lin_ϖ * x.τ') + @show size(dy.lin_τ * y.ϖ), size(dy.lin_ϖ * y.τ') + @show size(lin_τ .* ϖ') + lin_ϖ = ((dx.lin_τ * x.ϖ + + dx.lin_ϖ * x.τ') .+ + (dy.lin_τ * y.ϖ + + dy.lin_ϖ * y.τ') .- + lin_τ .* ϖ')./τ' + #end #@show xFT, xFT2, xFT3 - all(wx .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, y.Z⁺⁺, y.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dy.lin_Z⁺⁺, dy.lin_Z⁻⁺)) : nothing - all(wy .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dx.lin_Z⁺⁺, dx.lin_Z⁻⁺)) : nothing + #all(wx .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, y.Z⁺⁺, y.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dy.lin_Z⁺⁺, dy.lin_Z⁻⁺)) : nothing + #all(wy .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dx.lin_Z⁺⁺, dx.lin_Z⁻⁺)) : nothing - n = length(w); + n = length(τ); - wy = wy ./ w - wx = wx ./ w wx = reshape(wx,1,1,n) - wy = reshape(wy,1,1,n) + #wy = reshape(wy,1,1,n) - Z⁺⁺ = (wx .* xZ⁺⁺ .+ wy .* yZ⁺⁺) - Z⁻⁺ = (wx .* xZ⁻⁺ .+ wy .* yZ⁻⁺) - - tmpx1 = (dx.lin_τ.*x.ϖ .+ x.τ.*dx.lin_ϖ)./w - tmpy1 = (dy.lin_τ.*y.ϖ .+ y.τ.*dy.lin_ϖ)./w - tmpx2 = x.τ.*x.ϖ./w - tmpy2 = y.τ.*y.ϖ./w - xlin_Z⁺⁺ = tmpx1.*(x.Z⁺⁺.-Z⁺⁺) .+ tmpx2.*dx.dZ⁺⁺ - ylin_Z⁺⁺ = tmpy1.*(y.Z⁺⁺.-Z⁺⁺) .+ tmpy2.*dy.dZ⁺⁺ - xlin_Z⁻⁺ = tmpx1.*(x.Z⁻⁺.-Z⁻⁺) .+ tmpx2.*dx.dZ⁻⁺ - ylin_Z⁻⁺ = tmpy1.*(y.Z⁻⁺.-Z⁻⁺) .+ tmpy2.*dy.dZ⁻⁺ - lin_Z⁺⁺ = [xlin_Z⁺⁺, ylin_Z⁺⁺] - lin_Z⁻⁺ = [xlin_Z⁻⁺, ylin_Z⁻⁺] - - CoreScatteringOpticalProperties(τ, ϖ, Z⁺⁺, Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺) + Z⁺⁺ = (xZ⁺⁺ .* wx .+ wy .* yZ⁺⁺)./(wx .+ wy) + Z⁻⁺ = (xZ⁻⁺ .* wx .+ wy .* yZ⁻⁺)./(wx .+ wy) + + tmpx1 = (dx.lin_τ .* x.ϖ .+ dx.lin_ϖ .* x.τ') + tmpy1 = (dy.lin_τ .* y.ϖ .+ dy.lin_ϖ .* y.τ') + tmpx2 = reshape(wx, 1,1,1,n)#x.τ * x.ϖ #wx + tmpy2 = wy #y.τ * y.ϖ #wy + @show size(Z⁺⁺) + lin_Z⁺⁺ = zeros(nparams,size(Z⁺⁺,1), size(Z⁺⁺,2), size(Z⁺⁺,3)) + @show '1' + #ylin_Z⁺⁺ = zeros(nparams,size(Z⁺⁺,1), size(Z⁺⁺,2), size(Z⁺⁺,3)) + #@show '2' + lin_Z⁻⁺ = zeros(nparams,size(Z⁺⁺,1), size(Z⁺⁺,2), size(Z⁺⁺,3)) + @show '3' + #ylin_Z⁻⁺ = zeros(nparams,size(Z⁺⁺,1), size(Z⁺⁺,2), size(Z⁺⁺,3)) + #@show '4' + for i = 1:nparams + lin_Z⁺⁺[i,:,:,:] = (x.Z⁺⁺ .- Z⁺⁺).*reshape(tmpx1[i,:],1,1,n) .+ + dx.dZ⁺⁺[i,:,:].*tmpx2 .+ + y.Z⁺⁺.*reshape(tmpy1[i,:],1,1,n) .+ + dy.dZ⁺⁺[i,:,:].*tmpy2 + lin_Z⁻⁺[i,:,:,:] = (x.Z⁻⁺ .- Z⁻⁺).*reshape(tmpx1[i,:],1,1,n) .+ + dx.dZ⁻⁺[i,:,:].*tmpx2 .+ + y.Z⁻⁺.*reshape(tmpx1[i,:],1,1,n) .+ + dy.dZ⁻⁺[i,:,:].*tmpy2 + end + lin_Z⁺⁺ = lin_Z⁺⁺./reshape(wx .+ wy, 1,1,1,n) + lin_Z⁻⁺ = lin_Z⁻⁺./reshape(wx .+ wy, 1,1,1,n) + + return scatt_paired_xdx(CoreScatteringOpticalProperties(τ, ϖ, Z⁺⁺, Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺)) end # Concatenate Core Optical Properties, can have mixed dimensions! -function Base.:*( x::CoreScatteringOpticalProperties, y::CoreScatteringOpticalProperties, dx::linCoreScatteringOpticalProperties, dy::linCoreScatteringOpticalProperties ) +# Check with Christian what exactly this is supposed to do +#= +function Base.:*(x::CoreScatteringOpticalProperties, y::CoreScatteringOpticalProperties, dx::linCoreScatteringOpticalProperties, dy::linCoreScatteringOpticalProperties) arr_type = array_type(architecture(x.τ)) x = expandOpticalProperties(x,arr_type); y = expandOpticalProperties(y,arr_type); @@ -733,55 +800,51 @@ function Base.:*( x::CoreScatteringOpticalProperties, y::CoreScatteringOpticalPr CoreScatteringOpticalProperties([x.τ; y.τ],[x.ϖ; y.ϖ],cat(x.Z⁺⁺,y.Z⁺⁺, dims=3), cat(x.Z⁻⁺,y.Z⁻⁺, dims=3) ), linCoreScatteringOpticalProperties(cat(dx.dτ, dy.dτ, dims=2),cat(dx.dϖ, dy.dϖ, dims=2), cat(dx.dZ⁺⁺,dy.dZ⁺⁺, dims=4), cat(dx.dZ⁻⁺,dy.dZ⁻⁺, dims=4) ) end +=# + +function Base.:+(z1::scatt_paired_xdx, z2::abs_paired_xdx) + #(x::CoreScatteringOpticalProperties, dx::linCoreScatteringOpticalProperties, y::CoreAbsorptionOpticalProperties, dy::linCoreAbsorptionOpticalProperties) + # Predefine some arrays: + x = z1.x + dx = z1.dx + y = z2.x + dy = z2.dx -function Base.:+( x::CoreScatteringOpticalProperties, y::CoreAbsorptionOpticalProperties, dx::linCoreScatteringOpticalProperties, dy::linCoreAbsorptionOpticalProperties ) τ = x.τ .+ y.τ - wx = x.τ .* x.ϖ - ϖ = (wx) ./ τ + wx = x.τ .* x.ϖ + ϖ = wx ./ τ - lin_τ = [dx.lin_τ, dy.lin_τ] + lin_τ = dx.lin_τ + dy.lin_τ + lin_ϖ = (dx.lin_τ .* x.ϖ + x.τ .* dx.lin_ϖ - + ϖ .* lin_τ)./τ - lin_wx = (dx.lin_τ .* (x.ϖ - ϖ) + x.τ .* dx.lin_ϖ)./τ - lin_wy = -(dy.lin_τ .* ϖ)./τ - lin_ϖ = [lin_wx, lin_wy] - - tmpx1 = (dx.lin_τ.*x.ϖ .+ x.τ.*dx.lin_ϖ)./w - tmpx2 = x.τ.*x.ϖ./w - xlin_Z⁺⁺ = tmpx1.*(x.Z⁺⁺.-Z⁺⁺) .+ tmpx2.*dx.dZ⁺⁺ - xlin_Z⁻⁺ = tmpx1.*(x.Z⁻⁺.-Z⁻⁺) .+ tmpx2.*dx.dZ⁻⁺ - - lin_Z⁺⁺ = [xlin_Z⁺⁺, zeros(shape(Z⁺⁺))] #Check for size - lin_Z⁻⁺ = [xlin_Z⁻⁺, zeros(shape(Z⁻⁺))] + #@show xFT, xFT2, xFT3 + #all(wx .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, y.Z⁺⁺, y.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dy.lin_Z⁺⁺, dy.lin_Z⁻⁺)) : nothing + #all(wy .== 0.0) ? (return CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dx.lin_Z⁺⁺, dx.lin_Z⁻⁺)) : nothing - CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺) + return scatt_paired_xdx(CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dx.lin_Z⁺⁺, dx.lin_Z⁻⁺)) end -function Base.:+( y::CoreAbsorptionOpticalProperties, x::CoreScatteringOpticalProperties, dy::linCoreAbsorptionOpticalProperties, dx::linCoreScatteringOpticalProperties ) - τ = x.τ .+ y.τ - wx = x.τ .* x.ϖ - ϖ = (wx) ./ τ - - lin_τ = [dy.lin_τ, dx.lin_τ] - - lin_wx = (dx.lin_τ .* (x.ϖ - ϖ) + x.τ .* dx.lin_ϖ)./τ - lin_wy = -(dy.lin_τ .* ϖ)./τ - lin_ϖ = [lin_wy, lin_wx] - - tmpx1 = (dx.lin_τ.*x.ϖ .+ x.τ.*dx.lin_ϖ)./w - tmpx2 = x.τ.*x.ϖ./w +function Base.:+(z2::abs_paired_xdx, z1::scatt_paired_xdx) #(y::CoreAbsorptionOpticalProperties, dy::linCoreAbsorptionOpticalProperties, x::CoreScatteringOpticalProperties, dx::linCoreScatteringOpticalProperties) + x = z1.x + dx = z1.dx + y = z2.x + dy = z2.dx - xlin_Z⁺⁺ = tmpx1.*(x.Z⁺⁺.-Z⁺⁺) .+ tmpx2.*dx.dZ⁺⁺ - xlin_Z⁻⁺ = tmpx1.*(x.Z⁻⁺.-Z⁻⁺) .+ tmpx2.*dx.dZ⁻⁺ + τ = x.τ .+ y.τ + wx = x.τ .* x.ϖ + ϖ = wx ./ τ - lin_Z⁺⁺ = [zeros(shape(Z⁺⁺)), xlin_Z⁺⁺] #Check for size - lin_Z⁻⁺ = [zeros(shape(Z⁺⁺)), xlin_Z⁻⁺] + lin_τ = dx.lin_τ + dy.lin_τ + lin_ϖ = (dx.lin_τ .* x.ϖ + x.τ .* dx.lin_ϖ - + ϖ .* lin_τ)./τ - CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, lin_Z⁺⁺, lin_Z⁻⁺) + return scatt_paired_xdx(CoreScatteringOpticalProperties(τ, ϖ, x.Z⁺⁺, x.Z⁻⁺), linCoreScatteringOpticalProperties(lin_τ, lin_ϖ, dx.lin_Z⁺⁺, dx.lin_Z⁻⁺)) end - +#= function Base.:*( x::FT, y::CoreScatteringOpticalProperties{FT}, dy::linCoreScatteringOpticalProperties{FT} ) where FT CoreScatteringOpticalProperties(y.τ * x, y.ϖ, y.Z⁺⁺, y.Z⁻⁺, y.G), linCoreScatteringOpticalProperties(dy.lin_τ * x, dy.lin_ϖ, dy.linZ⁺⁺, dy.linZ⁻⁺, dy.linG) end - +=# diff --git a/src/CoreRT/parameters_from_yaml.jl b/src/CoreRT/parameters_from_yaml.jl index 3b2e31eb..730a04cc 100644 --- a/src/CoreRT/parameters_from_yaml.jl +++ b/src/CoreRT/parameters_from_yaml.jl @@ -7,7 +7,11 @@ the input, and producing the output object. =# "Check that a field exists in yaml file" -function check_yaml_field(dict::Dict{Any, Any}, full_keys::Array{String}, curr_keys::Array{String}, final_type::Type, valid_options::Array{String}) +function check_yaml_field(dict::Dict{Any, Any}, + full_keys::Array{String}, + curr_keys::Array{String}, + final_type::Type, + valid_options::Array{String}) # Should have at least one key @assert length(curr_keys) >= 1 @@ -38,8 +42,8 @@ function validate_aerosols(aerosols::Union{Array{Dict{Any, Any}}, Vector{Any}}) (["σ"], Real), (["nᵣ"], Real), (["nᵢ"], Real), - (["p₀"], Real), - (["p₀"], Real)] + (["z₀"], Real), # (["p₀"], Real), (Suniti) + (["σ₀"], Real)] # (["p₀"], Real)] (Suniti) for aerosol in aerosols for i in 1:length(fields) @@ -62,8 +66,8 @@ function aerosol_params_to_obj(aerosols::Union{Array{Dict{Any, Any}}, Vector{Any FT(aerosol["nᵣ"]), FT(aerosol["nᵢ"])) - new_rt_aerosol_obj = RT_Aerosol(new_aerosol_obj, FT(aerosol["τ_ref"]), Normal(FT(aerosol["p₀"]), FT(aerosol["σp"]))) - + # new_rt_aerosol_obj = RT_Aerosol(new_aerosol_obj, FT(aerosol["τ_ref"]), Normal(FT(aerosol["p₀"]), FT(aerosol["σp"]))) + new_rt_aerosol_obj = RT_Aerosol(new_aerosol_obj, FT(aerosol["τ_ref"]), LogNormal(FT(aerosol["z₀"]), FT(aerosol["σ₀"]))) # Suniti push!(rt_aerosol_obj_list, new_rt_aerosol_obj) end diff --git a/src/CoreRT/rt_helper_functions.jl b/src/CoreRT/rt_helper_functions.jl index f5040a75..75c94345 100644 --- a/src/CoreRT/rt_helper_functions.jl +++ b/src/CoreRT/rt_helper_functions.jl @@ -73,24 +73,44 @@ end default_matrix(FT, arr_type, dims, nSpec) = arr_type(zeros(FT, tuple(dims[1], dims[2], nSpec))) "Default matrix in ieRT calculation (zeros)" default_matrix_ie(FT, arr_type, dims, nSpec, nRaman) = arr_type(zeros(FT, tuple(dims[1], dims[2], nSpec, nRaman))) +"Default Jacobian matrix in RT calculation (zeros)" +default_matrix_lin(FT, arr_type, nparams, dims, nSpec) = arr_type(zeros(FT, tuple(nparams, dims[1], dims[2], nSpec))) +"Default Jacobian matrix in ieRT calculation (zeros)" +default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, nRaman) = arr_type(zeros(FT, tuple(nparams, dims[1], dims[2], nSpec, nRaman))) "Default J matrix in RT calculation (zeros)" default_J_matrix(FT, arr_type, dims, nSpec) = arr_type(zeros(FT, tuple(dims[1], 1, nSpec))) "Default J matrix in ieRT calculation (zeros)" default_J_matrix_ie(FT, arr_type, dims, nSpec, nRaman) = arr_type(zeros(FT, tuple(dims[1], 1, nSpec, nRaman))) +"Default J' matrix in RT calculation (zeros)" +default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec) = arr_type(zeros(FT, tuple(nparams, dims[1], 1, nSpec))) +"Default J' matrix in ieRT calculation (zeros)" +default_J_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, nRaman) = arr_type(zeros(FT, tuple(nparams, dims[1], 1, nSpec, nRaman))) "Default matrix in RT calculation (zeros)" default_matrix(FT, arr_type, NSens, dims, nSpec) = [arr_type(zeros(FT, (dims[1], dims[2], nSpec))) for i=1:NSens] #arr_type(zeros(FT, tuple(NSens, dims[1], dims[2], nSpec))) "Default matrix in ieRT calculation (zeros)" -default_matrix_ie(FT, arr_type, NSens, dims, nSpec, nRaman) = [zeros(FT, (dims[1], dims[2], nSpec, nRaman)) for i=1:NSens] +default_matrix_ie(FT, arr_type, NSens, dims, nSpec, nRaman) = [arr_type(zeros(FT, (dims[1], dims[2], nSpec, nRaman))) for i=1:NSens] +#zeros(FT, tuple(NSens, dims[1], dims[2], nSpec, nRaman))) +"Default Jacobian matrix in RT calculation (zeros)" +default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec) = [arr_type(zeros(FT, (nparams, dims[1], dims[2], nSpec))) for i=1:NSens] +#arr_type(zeros(FT, tuple(NSens, dims[1], dims[2], nSpec))) +"Default Jacobian matrix in ieRT calculation (zeros)" +default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, nRaman) = [arr_type(zeros(FT, (nparams, dims[1], dims[2], nSpec, nRaman))) for i=1:NSens] #zeros(FT, tuple(NSens, dims[1], dims[2], nSpec, nRaman))) "Default J matrix in RT calculation (zeros)" default_J_matrix(FT, arr_type, NSens, dims, nSpec) = [arr_type(zeros(FT, (dims[1], 1, nSpec))) for i=1:NSens] #arr_type(zeros(FT, tuple(NSens, dims[1], 1, nSpec))) "Default J matrix in ieRT calculation (zeros)" -default_J_matrix_ie(FT, arr_type, NSens, dims, nSpec, nRaman) = [zeros(FT, (dims[1], 1, nSpec, nRaman)) for i=1:NSens] +default_J_matrix_ie(FT, arr_type, NSens, dims, nSpec, nRaman) = [arr_type(zeros(FT, (dims[1], 1, nSpec, nRaman))) for i=1:NSens] +#arr_type(zeros(FT, tuple(NSens, dims[1], 1, nSpec, nRaman))) +"Default J' matrix in RT calculation (zeros)" +default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec) = [arr_type(zeros(FT, (nparams, dims[1], 1, nSpec))) for i=1:NSens] +#arr_type(zeros(FT, tuple(NSens, dims[1], 1, nSpec))) +"Default J' matrix in ieRT calculation (zeros)" +default_J_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, nRaman) = [arr_type(zeros(FT, (nparams, dims[1], 1, nSpec, nRaman))) for i=1:NSens] #arr_type(zeros(FT, tuple(NSens, dims[1], 1, nSpec, nRaman))) ##### Only for testing, random matrices: @@ -219,51 +239,54 @@ make_lin_added_layer(RS_type::Union{noRS, noRS_plus}, nSpec, nparams, redn1, redn2) = linAddedLayer( - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, redn1, dims, nSpec), - default_matrix(FT, arr_type, redn1, dims, nSpec), - default_matrix(FT, arr_type, redn1, dims, nSpec), - default_matrix(FT, arr_type, redn1, dims, nSpec), - default_J_matrix(FT, arr_type, redn2, dims, nSpec), - default_J_matrix(FT, arr_type, redn2, dims, nSpec) + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_J_matrix_lin(FT, arr_type, redn2, dims, nSpec), + default_J_matrix_lin(FT, arr_type, redn2, dims, nSpec) ) "Make an added layer, supplying all default matrices" -make_lin_added_layer(RS_type::Union{RRS, RRS_plus,VS_0to1_plus, VS_1to0_plus}, lin::vSmartMOM_lin, FT, arr_type, dims, nSpec, nparams) = linAddedLayerRS( - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix(FT, arr_type, 3, dims, nSpec), - default_matrix(FT, arr_type, 3, dims, nSpec), - default_matrix(FT, arr_type, 3, dims, nSpec), - default_matrix(FT, arr_type, 3, dims, nSpec), - default_J_matrix(FT, arr_type, 4, dims, nSpec), - default_J_matrix(FT, arr_type, 4, dims, nSpec), - default_matrix_ie(FT, arr_type, 3, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, 3, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, 3, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, 3, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, 4, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, 4, dims, nSpec, RS_type.n_Raman) - ) +make_lin_added_layer(RS_type::Union{RRS, RRS_plus,VS_0to1_plus, VS_1to0_plus}, + lin::vSmartMOM_lin, + FT, arr_type, dims, nSpec, nparams, redn1, redn2) = + linAddedLayerRS( + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_matrix_lin(FT, arr_type, redn1, dims, nSpec), + default_J_matrix_lin(FT, arr_type, redn2, dims, nSpec), + default_J_matrix_lin(FT, arr_type, redn2, dims, nSpec), + default_matrix_ie_lin(FT, arr_type, redn1, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, redn1, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, redn1, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, redn1, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, redn2, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, redn2, dims, nSpec, RS_type.n_Raman) + ) - +#= "Make a random added layer, supplying all random matrices" -make_lin_added_layer_rand(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, FT, arr_type, dims, nSpec, naparams) = linAddedLayer( +make_lin_added_layer_rand(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, FT, arr_type, dims, nSpec, nparams) = linAddedLayer( default_matrix_rand(FT, arr_type, nparams, dims, nSpec), default_matrix_rand(FT, arr_type, nparams, dims, nSpec), default_matrix_rand(FT, arr_type, nparams, dims, nSpec), @@ -277,7 +300,7 @@ make_lin_added_layer_rand(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, F default_J_matrix_rand(FT, arr_type, 4, dims, nSpec), default_J_matrix_rand(FT, arr_type, 4, dims, nSpec) ) - +=# "Make a composite layer, supplying all default matrices" make_lin_composite_layer(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, FT, arr_type, dims, nSpec, nparams) = linCompositeLayer( @@ -287,27 +310,27 @@ make_lin_composite_layer(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, #default_matrix(FT, arr_type, dims, nSpec), #default_J_matrix(FT, arr_type, dims, nSpec), #default_J_matrix(FT, arr_type, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec)) + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec)) "Make a composite layer, supplying all default matrices" make_lin_composite_layer(RS_type::Union{RRS, RRS_plus,VS_0to1_plus, VS_1to0_plus}, lin::vSmartMOM_lin, FT, arr_type, dims, nSpec, nparams) = linCompositeLayerRS( - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, dims, nSpec), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman) + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, dims, nSpec), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, dims, nSpec, RS_type.n_Raman) ) "Make a composite layer, supplying all default matrices" @@ -326,18 +349,18 @@ make_lin_composite_layer(RS_type::Union{noRS, noRS_plus}, lin::vSmartMOM_lin, #default_matrix(FT, arr_type, NSens, dims, nSpec), #default_J_matrix(FT, arr_type, NSens, dims, nSpec), #default_J_matrix(FT, arr_type, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec)) + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec)) "Make a composite layer, supplying all default matrices" make_lin_composite_layer(RS_type::Union{RRS, RRS_plus, VS_0to1_plus, VS_1to0_plus}, lin::vSmartMOM_lin, @@ -367,30 +390,30 @@ make_lin_composite_layer(RS_type::Union{RRS, RRS_plus, VS_0to1_plus, VS_1to0_plu #default_matrix_ie(FT, arr_type, NSens, dims, nSpec, RS_type.n_Raman), #default_J_matrix_ie(FT, arr_type, NSens, dims, nSpec, RS_type.n_Raman), #default_J_matrix_ie(FT, arr_type, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_J_matrix(FT, arr_type, nparams, NSens, dims, nSpec), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), - default_J_matrix_ie(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman) + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_J_matrix_lin(FT, arr_type, nparams, NSens, dims, nSpec), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman), + default_J_matrix_ie_lin(FT, arr_type, nparams, NSens, dims, nSpec, RS_type.n_Raman) ) #ending linearization block diff --git a/src/CoreRT/rt_run.jl b/src/CoreRT/rt_run.jl index db2f7db1..14a4fcbb 100644 --- a/src/CoreRT/rt_run.jl +++ b/src/CoreRT/rt_run.jl @@ -209,13 +209,17 @@ function rt_run(RS_type::AbstractRamanType, end function rt_run_test(RS_type::AbstractRamanType, - model::vSmartMOM_Model, lin::vSmartMOM_lin, nparams::Int, + model::vSmartMOM_Model, + lin::vSmartMOM_lin, + nparams::Int, iBand) rt_run(RS_type, model, lin, nparams, iBand) end function rt_run(RS_type::AbstractRamanType, - model::vSmartMOM_Model, lin::vSmartMOM_lin, nparams::Int, iBand) + model::vSmartMOM_Model, + lin::vSmartMOM_lin, + nparams::Int, iBand) @unpack obs_alt, sza, vza, vaz = model.obs_geom # Observational geometry properties @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, μ₀, iμ₀, Nquad = model.quad_points # All quadrature points pol_type = model.params.polarization_type @@ -261,7 +265,7 @@ function rt_run(RS_type::AbstractRamanType, # β_trunc, and Z. After that, however, it is # necessary to account for all parameters # explicitly using the chain rule to account - # for the varying dependence of τ, ϖ, β_trunc, + # for the varying dependence of τ, ϖ, β_trunc, # and Z on the parameters for different # atmospheric layers. Note: in cases of weak # atmospheric absorption, this problem can be @@ -275,8 +279,8 @@ function rt_run(RS_type::AbstractRamanType, # 4. nᵢ x Naer # 5. r₀ x Naer # 6. σ₀ x Naer - # 7. Pᵥ x Naer - # 8. σᵥ x Naer + # 7. z₀ x Naer + # 8. σz x Naer # 9. Psurf # 10. VMR1_CO2 (for f₁(P)) # 11. VMR2_CO2 (for f₂(P)) @@ -316,6 +320,7 @@ function rt_run(RS_type::AbstractRamanType, # Create arrays @timeit "Creating layers" added_layer = make_added_layer(RS_type, FT_dual, arr_type, dims, nSpec) + @show dims, nSpec, nparams @timeit "Creating layers" lin_added_layer = make_lin_added_layer(RS_type, lin, FT, arr_type, dims, nSpec, nparams, Int64(3), Int64(4)) @@ -351,13 +356,13 @@ function rt_run(RS_type::AbstractRamanType, InelasticScattering.computeRamanZλ!(RS_type, pol_type,Array(qp_μ), m, arr_type) # Compute the core layer optical properties: @timeit "OpticalProps" layer_opt_props, fScattRayleigh, - lin_layer_opt_props, lin_fScattRayleigh = + lin_layer_opt_props = constructCoreOpticalProperties(RS_type, iBand, m, model, lin); # Determine the scattering interface definitions: scattering_interfaces_all, τ_sum_all, lin_τ_sum_all = extractEffectiveProps(layer_opt_props, lin_layer_opt_props, quad_points); #@show typeof(layer_opt_props) - nparams = length(lin_fScattRayleigh[:,1]) + #nparams = length(lin_fScattRayleigh[:,1]) speclen = length(RS_type.fscattRayl[1,:]) # Loop over vertical layers: @showprogress 1 "Looping over layers ..." for iz = 1:Nz # Count from TOA to BOA @@ -368,19 +373,18 @@ function rt_run(RS_type::AbstractRamanType, if !(typeof(RS_type) <: noRS) RS_type.fscattRayl = expandBandScalars(RS_type, fScattRayleigh[iz]) - - RS_type.lin_fscattRayl = zeros(nparams,speclen) - for ctr = 1:nparams - RS_type.lin_fScattRayleigh[ctr,:] = - expandBandScalars(RS_type, lin_fScattRayleigh[ctr,iz]) - end + #RS_type.lin_fscattRayl = zeros(nparams,speclen) + #for ctr = 1:nparams + # RS_type.lin_fScattRayleigh[ctr,:] = + # expandBandScalars(RS_type, lin_fScattRayleigh[ctr,iz]) + #end end # Expand all layer optical properties to their full dimension: @timeit "OpticalProps" layer_opt = expandOpticalProperties(layer_opt_props[iz], arr_type) @timeit "OpticalProps" lin_layer_opt = - expandOpticalProperties(lin_layer_opt_props[iz], arr_type) + expandOpticalProperties_lin(lin_layer_opt_props[iz], arr_type) # Perform Core RT (doubling/elemental/interaction) rt_kernel!(RS_type, pol_type, SFI, @@ -390,8 +394,8 @@ function rt_run(RS_type::AbstractRamanType, layer_opt, lin_layer_opt, scattering_interfaces_all[iz], - τ_sum_all[:,iz], - lin_τ_sum_all[:,:,iz], + τ_sum_all[:,iz], μ₀, + #lin_τ_sum_all[:,:,iz], m, quad_points, I_static, model.params.architecture, diff --git a/src/Scattering/Scattering.jl b/src/Scattering/Scattering.jl index 092dde95..a4673551 100644 --- a/src/Scattering/Scattering.jl +++ b/src/Scattering/Scattering.jl @@ -55,5 +55,6 @@ export compute_B, compute_ab, GreekCoefs, comp_ab, compute_mie_π_τ!, ConjugateTransposePairs, AbstractPolarizationType, AbstractAerosolType, AbstractAerosolType, MieModel, AbstractTruncationType, phase_function +export compute_ref_aerosol_extinction_lin, compute_aerosol_optical_properties_lin end diff --git a/src/Scattering/lin_compute_NAI2.jl b/src/Scattering/lin_compute_NAI2.jl index 6802b5de..292195e6 100644 --- a/src/Scattering/lin_compute_NAI2.jl +++ b/src/Scattering/lin_compute_NAI2.jl @@ -14,7 +14,7 @@ Input: MieModel, holding all computation and aerosol properties Output: AerosolOptics, holding all Greek coefficients and Cross-Sectional information """ #TODO: define linAerosolOptics -function compute_aerosol_optical_properties(model::MieModel{FDT}, FT2::Type=Float64) where FDT <: NAI2 +function compute_aerosol_optical_properties_lin(model::MieModel{FDT}, FT2::Type=Float64) where FDT <: NAI2 # Unpack the model @unpack aerosol, λ, r_max, nquad_radius = model @@ -324,23 +324,23 @@ function compute_aerosol_optical_properties(model::MieModel{FDT}, FT2::Type=Floa dAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=dω̃, dk=d_bulk_C_ext, dk_ref=zeros(FT, 4), dfᵗ=zeros(FT, 4)) else =# - greek_coefs = GreekCoefs(α,β,γ,δ,ϵ,ζ) - #d_greek_coefs=[]; - #for i=1:4 - d_greek_coefs = [GreekCoefs(dα[i,:], - dβ[i,:], - dγ[i,:], - dδ[i,:], - dϵ[i,:], - dζ[i,:]) for i=1:4] - #push!(d_greek_coefs, tmp_greek_coefs); - #end - return AerosolOptics(greek_coefs=greek_coefs, ω̃=ω̃, k=(bulk_C_ext), k_ref=FT(0), fᵗ=FT(1)), - dAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=dω̃, dk=d_bulk_C_ext, dk_ref=zeros(FT, 4), dfᵗ=zeros(FT, 4)) + greek_coefs = GreekCoefs(α,β,γ,δ,ϵ,ζ) + #d_greek_coefs=[]; + #for i=1:4 + d_greek_coefs = [GreekCoefs(dα[i,:], + dβ[i,:], + dγ[i,:], + dδ[i,:], + dϵ[i,:], + dζ[i,:]) for i=1:4] + #push!(d_greek_coefs, tmp_greek_coefs); #end + return AerosolOptics(greek_coefs=greek_coefs, ω̃=ω̃, k=(bulk_C_ext), k_ref=FT(0), fᵗ=FT(1)), + dAerosolOptics(d_greek_coefs=d_greek_coefs, dω̃=dω̃, dk=d_bulk_C_ext, dk_ref=zeros(FT, 4), dfᵗ=zeros(FT, 4)) end -function compute_ref_aerosol_extinction(model::MieModel{FDT}, FT2::Type=Float64) where FDT <: NAI2 +# function compute_ref_aerosol_extinction_lin(model::MieModel{FDT}, FT2::Type=Float64) where FDT <: NAI2 +function compute_ref_aerosol_extinction_lin(model::MieModel)#{FT}, FT) #2::Type=Float64) where FDT <: NAI2 # Unpack the model @unpack computation_type, aerosol, λ, polarization_type, r_max, nquad_radius = model @@ -440,7 +440,7 @@ end Compute phase function from aerosol distribution with log-normal mean μ [µm] and σ Output: μ, w_μ, P, C_ext, C_sca, g """ -function phase_function(aerosol::Aerosol, λ, r_max, nquad_radius) +function phase_function_lin(aerosol::Aerosol, λ, r_max, nquad_radius) # Extract variables from aerosol struct: @unpack size_distribution, nᵣ, nᵢ = aerosol @@ -610,7 +610,7 @@ end Compute phase function from mono-modal aerosol with radius `r` at wavelength `λ`, both in `[μm]` Output: μ, w_μ, P, C_ext, C_sca, g """ -function phase_function(r::FT, λ::FT, nᵣ::FT, nᵢ::FT) where {FT<:AbstractFloat} +function phase_function_lin(r::FT, λ::FT, nᵣ::FT, nᵢ::FT) where {FT<:AbstractFloat} # Imaginary part of the refractive index must be ≥ 0 (definition) @assert nᵢ ≥ 0 "Imaginary part of the refractive index must be ≥ 0 (definition)" # Wavenumber diff --git a/src/Scattering/lin_mie_helper_functions.jl b/src/Scattering/lin_mie_helper_functions.jl index 0e166594..7e1f441c 100644 --- a/src/Scattering/lin_mie_helper_functions.jl +++ b/src/Scattering/lin_mie_helper_functions.jl @@ -90,7 +90,7 @@ Compute all an, bn using compute_mie_ab! Input: MieModel, wavelength (λ), radius Output: an, bn. Both of shape (aerosol.nquad_radius, N_max) (N_max from aerosol.r_max) """ -function compute_anbn(model::MieModel, λ, radius) +function compute_anbn_lin(model::MieModel, λ, radius) @unpack computation_type, aerosol, r_max, nquad_radius, λ, polarization_type, truncation_type, wigner_A, wigner_B = model @unpack size_distribution, nᵣ, nᵢ = aerosol @@ -136,12 +136,12 @@ end From the an, bn matrices, precompute all (an✶)am, (an✶)bm, (bn✶)am, (bn✶)bm This allows quick computation of (an✶ + bn✶) × (am + bm) """ -function compute_avg_anbns!(an, bn, ab_pairs, w, Nmax, N_max_) +function compute_avg_anbns!(an, bn, ab_pairs, dan, dbn, d_ab_pairs, w, Nmax, N_max_) FT2 = eltype(an) # Unpack ab_pairs mat_anam, mat_anbm, mat_bnam, mat_bnbm = ab_pairs - mat_danam, mat_danbm, mat_dbnam, mat_dbnbm = ab_pairs + mat_danam, mat_danbm, mat_dbnam, mat_dbnbm = d_ab_pairs # Fill all matrices with 0 [fill!(mat, 0) for mat in [mat_anam, mat_bnbm, mat_anbm, mat_bnam]] [fill!(mat, 0) for mat in [mat_danam, mat_dbnbm, mat_danbm, mat_dbnam]] @@ -271,7 +271,7 @@ function reconstruct_phase(greek_coefs, dgreek_coefs, μ; returnLeg=false) # For truncation in δ-BGE, we need P and P² as well, convenient to return here: return returnLeg ? (scattering_matrix, dscattering_matrix, P, P²) : (scattering_matrix, dscattering_matrix) end - +#= """ $(FUNCTIONNAME)(depol) Returns the greek coefficients (as [`GreekCoefs`](@ref)) of the Rayleigh phase function given @@ -293,7 +293,7 @@ function get_greek_rayleigh(depol::Number) ζ = FT[0.0, 0.0, 0.0] return GreekCoefs(α, β, γ, δ, ϵ, ζ) end - +=# """ $(FUNCTIONNAME)(k, an, bn, w) Calculate the average Scattering and Extinction Cross Section @@ -303,10 +303,13 @@ function compute_avg_C_scatt_ext(k, an, bn, dan, dbn, w) n_ = collect(1:size(an)[2]); n_ = 2n_ .+ 1 coef = 2π / k^2 * n_' - return (coef * (w' * (abs2.(an') + abs2.(bn'))')', coef * (w' * real(an + bn))', - 2*coef * (w' * (real(an)*real(dan') + real(bn)*real(dbn') + imag(an)*imag(dan') + imag(bn)*imag(dbn'))')', coef * (w' * real(dan + dbn))') + return (coef * (w' * (abs2.(an') + abs2.(bn'))')', + coef * (w' * real(an + bn))', + 2*coef * (w' * (real(an)*real(dan') + real(bn)*real(dbn') + imag(an)*imag(dan') + imag(bn)*imag(dbn'))')', + coef * (w' * real(dan + dbn))') end - + +#= """ Compute probability weights of radii """ function compute_wₓ(size_distribution, wᵣ, r, r_max) @@ -393,7 +396,7 @@ See Sanghavi 2014, eq. 16 construct_B_matrix(mod::Stokes_I, α, β, γ, δ, ϵ, ζ, l::Int) = β[l] -#= + """ $(FUNCTIONNAME)(mod::AbstractPolarizationType, μ, α, β, γ, δ, ϵ, ζ, m::Int) Compute moments of the phase matrix diff --git a/src/Scattering/lin_types.jl b/src/Scattering/lin_types.jl index e224d16f..c3f3a40f 100644 --- a/src/Scattering/lin_types.jl +++ b/src/Scattering/lin_types.jl @@ -75,7 +75,7 @@ A struct which holds all computed aerosol optics $(DocStringExtensions.FIELDS) """ #Base.@kwdef -Base.@kwdef mutable struct dAerosolOptics{FT<:AbstractFloat} +Base.@kwdef mutable struct dAerosolOptics{FT}#<:AbstractFloat} "derivatives of Greek matrix w.r.t nᵣ, nᵢ, r₀ and σ₀" d_greek_coefs::Vector{GreekCoefs{FT}} "derivatives of Single Scattering Albedo w.r.t nᵣ, nᵢ, r₀ and σ₀" diff --git a/src/vSmartMOM.jl b/src/vSmartMOM.jl index 570ca52e..7c76cd41 100644 --- a/src/vSmartMOM.jl +++ b/src/vSmartMOM.jl @@ -47,7 +47,7 @@ include("SolarModel/SolarModel.jl") # Export some vSmartMOM functions export default_parameters, parameters_from_yaml, model_from_parameters, rt_run - +export lin_model_from_parameters using .Architectures using .Absorption diff --git a/test/prototyping/Balsamic_OCO2_test.jl b/test/prototyping/Balsamic_OCO2_test.jl index 0780f4a6..db1d00e6 100644 --- a/test/prototyping/Balsamic_OCO2_test.jl +++ b/test/prototyping/Balsamic_OCO2_test.jl @@ -214,6 +214,18 @@ end parameters = parameters_from_yaml("test/test_parameters/3BandParameters.yaml") #parameters.architecture = CPU() FT = Float64 +RS_type = InelasticScattering.noRS( + fscattRayl = [FT(1)], + ϖ_Cabannes = [FT(1)], + bandSpecLim = [], + iBand = [1], + F₀ = zeros(FT,1,1)); +#RS_type.F₀ = zeros(model.params.polarization_type.n, length(P)) +#for i=1:length(P) +# sol_trans = Tsolar_interp(ν[i]); +# F₀[i] = sol_trans * P[i]; +# RS_type.F₀[1,i] = F₀[i]; +#end # Load OCO Data: # File names: @@ -341,29 +353,42 @@ parameters.vaz = v_raz # Produce black-body in wavenumber range Tsolar = solar_transmission_from_file("/home/rjeyaram/vSmartMOM/src/SolarModel/solar.out") Tsolar_interp = LinearInterpolation(Tsolar[:, 1], Tsolar[:, 2]) - +# New parameter order: + # 1: p_surf + # (2-7) CO2: VMR₀, VMR₁, VMR₂, HCO2, ln(z₀CO2), σCO2 + # (8,9) H2O: a1, a2 + # (10-16) Aerosol1: ln(τ_ref), nᵣ, nᵢ, r₀, σᵣ, z₀, σz + # (17-23) Aerosol2: ln(τ_ref), nᵣ, nᵢ, r₀, σᵣ, z₀, σz + # (24-32) Albedo: ρA₁, ρA₂, ρA₃, ρW₁, ρW₂, ρW₃, ρS₁, ρS₂, ρS₃ + # BRDF: RPV/RossLi/CoxMunk +x = FT[mean(oco_soundings[1].p_half[end]), # 1: Psurf (hPa) - currently assumed to be the mean of all multiangle modes - think of a way to separate them eventually + 400.0e-6, 0.0, 0.0, 8.0, log(3.0), 0.62, # (2-7) CO2: VMR₀, VMR₁, VMR₂, HCO2, ln(z₀CO2), σCO2 + 0.0, 0.0, # (8,9) H2O: a1, a2 + -1.2, 1.42, 0.001, log(0.1), 0.83, log(5), 0.47, # (10-16) Aerosol1: ln(τ_ref), nᵣ, nᵢ, r₀, σᵣ, z₀, σz + -2.3, 1.65, 0.001, log(0.5), 0.60, log(3), 0.64, # (17-13) Aerosol1: ln(τ_ref), nᵣ, nᵢ, r₀, σᵣ, z₀, σz + 0.6, 0.0, 0.0, 0.4, 0.0, 0.0, 0.2, 0.0, 0.0] # (24-32) Albedo: ρA₁, ρA₂, ρA₃, ρW₁, ρW₂, ρW₃, ρS₁, ρS₂, ρS₃ #Initial guess state vector -x = FT[0.2377, # 1. Legendre Lambertian Albedo, A-band - -3.5, # log(AOD) - -0.00244, # 2. Legendre Lambertian Albedo, A-band - 0, # 3. Legendre Lambertian Albedo, A-band - 0, # 3. Legendre Lambertian Albedo, W-band - 0, # 3. Legendre Lambertian Albedo, S-band - 0.2, # 1. Legendre Lambertian Albedo, W-band - 0, # 2. Legendre Lambertian Albedo, W-band - 0.2, # 1. Legendre Lambertian Albedo, S-band - 0, # 2. Legendre Lambertian Albedo, S-band - 0, # H2O - 400e-6, # CO2 1 - 400e-6, # CO2 2 - 400e-6, # CO2 3 - 0, - 1.33, # nᵣ - 0.01, # nᵢ - -1.69, # lnr₀ - 0.3, # σ₀ - 690., #p₀ - 50.] #σₚ +#x = FT[0.2377, # 1. Legendre Lambertian Albedo, A-band +# -3.5, # log(AOD) +# -0.00244, # 2. Legendre Lambertian Albedo, A-band +# 0, # 3. Legendre Lambertian Albedo, A-band +# 0, # 3. Legendre Lambertian Albedo, W-band +# 0, # 3. Legendre Lambertian Albedo, S-band +# 0.2, # 1. Legendre Lambertian Albedo, W-band +# 0, # 2. Legendre Lambertian Albedo, W-band +# 0.2, # 1. Legendre Lambertian Albedo, S-band +# 0, # 2. Legendre Lambertian Albedo, S-band +# 0, # H2O +# 400e-6, # CO2 1 +# 400e-6, # CO2 2 +# 400e-6, # CO2 3 +# 0, +# 1.33, # nᵣ +# 0.01, # nᵢ +# -1.69, # lnr₀ +# 0.3, # σ₀ +# 690., #p₀ +# 50.] #σₚ nparams=length(x) # Run FW model: diff --git a/test/prototyping/runner.jl b/test/prototyping/runner.jl index 14f0695a..c468f47e 100644 --- a/test/prototyping/runner.jl +++ b/test/prototyping/runner.jl @@ -116,21 +116,42 @@ end function lin_runner!(y, x, parameters=parameters, oco_sounding= oco_soundings, Tsolar = Tsolar_interp) Nangles = length(oco_sounding); #@show Nangles + # New parameter order: + # 1: p_surf + # (2-7) CO2: VMR₀, VMR₁, VMR₂, HCO2, z₀CO2, σCO2 + # (8,9) H2O: a1, a2 + # (10-16) Aerosol1: ln(τ_ref), nᵣ, nᵢ, r₀, σᵣ, z₀, σz + # (17-23) Aerosol2: ln(τ_ref), nᵣ, nᵢ, r₀, σᵣ, z₀, σz + # (24-32) Albedo: ρA₁, ρA₂, ρA₃, ρW₁, ρW₂, ρW₃, ρS₁, ρS₂, ρS₃ + # BRDF: RPV/RossLi/CoxMunk # Set parameters fields as the dual numbers - parameters.brdf = [CoreRT.LambertianSurfaceLegendre([x[1],x[3],x[4]]), - CoreRT.LambertianSurfaceLegendre([x[7],x[8],x[5]]), - CoreRT.LambertianSurfaceLegendre([x[9],x[10],x[6]])]; + parameters.brdf = [CoreRT.LambertianSurfaceLegendre([x[24],x[25],x[26]]), + CoreRT.LambertianSurfaceLegendre([x[27],x[28],x[29]]), + CoreRT.LambertianSurfaceLegendre([x[30],x[31],x[32]])]; - parameters.scattering_params.rt_aerosols[1].τ_ref = exp(x[2]); - @show x[2], exp(x[2]) - parameters.scattering_params.rt_aerosols[1].profile = Normal(x[20], x[21]); #800.0; #x[4] - #parameters.scattering_params.rt_aerosols[1].profile.σp = x[21]; - @show x[20], x[21] - parameters.scattering_params.rt_aerosols[1].aerosol.size_distribution = LogNormal(x[18], x[19]); #x[4] + parameters.scattering_params.rt_aerosols[1].τ_ref = exp(x[10]); + @show x[10], exp(x[10]) + parameters.scattering_params.rt_aerosols[1].aerosol.nᵣ = x[11] + parameters.scattering_params.rt_aerosols[1].aerosol.nᵢ = x[12] + @show x[11], x[12] + parameters.scattering_params.rt_aerosols[1].aerosol.size_distribution = + LogNormal(x[13], x[14]); #x[4] + @show x[13], x[14] + parameters.scattering_params.rt_aerosols[1].profile = + LogNormal(x[15], x[16]); #800.0; #x[4] + @show x[15], x[16] + + parameters.scattering_params.rt_aerosols[2].τ_ref = exp(x[17]); + @show x[17], exp(x[17]) + parameters.scattering_params.rt_aerosols[2].aerosol.nᵣ = x[18] + parameters.scattering_params.rt_aerosols[2].aerosol.nᵢ = x[19] @show x[18], x[19] - parameters.scattering_params.rt_aerosols[1].aerosol.nᵣ = x[16] - parameters.scattering_params.rt_aerosols[1].aerosol.nᵢ = x[17] - @show x[16], x[17] + parameters.scattering_params.rt_aerosols[2].aerosol.size_distribution = + LogNormal(x[20], x[21]); #x[4] + @show x[20], x[21] + parameters.scattering_params.rt_aerosols[2].profile = + LogNormal(x[22], x[23]); #800.0; #x[4] + @show x[22], x[23] nparams=Int(length(x)) #parameters.p = oco_sounding.p_half @@ -138,12 +159,18 @@ function lin_runner!(y, x, parameters=parameters, oco_sounding= oco_soundings, T #parameters.T = oco_sounding.T# .+ 1.0 #.+ x[15] #parameters.sza = oco_sounding.sza #parameters.vza = [oco_sounding.vza] + psurf = x[1] for ia=1:Nangles #if ia==1 parameters.p = oco_sounding[ia].p_half + #inducing psurf + parameters.p[end] = psurf + #@show parameters.q#[end], parameters.p[end-1] + @assert( parameters.p[end] >= parameters.p[end-1] ) parameters.q = oco_sounding[ia].q parameters.T = oco_sounding[ia].T# .+ 1.0 #.+ x[15] parameters.sza = oco_sounding[ia].sza + #parameters.vza = [oco_sounding[ia].vza] #parameters.vaz = [oco_sounding[ia].raa] #else @@ -151,13 +178,36 @@ function lin_runner!(y, x, parameters=parameters, oco_sounding= oco_soundings, T # parameters.vaz = vcat(parameters.vaz, oco_sounding[ia].raa) #end end - parameters.absorption_params.vmr["H2O"] = [parameters.q[1:65]*x[11] * 1.8; - parameters.q[66:end]*x[15] * 1.8]; - a1 = zeros(7) .+ x[12] - a2 = zeros(7) .+ x[13] - a3 = zeros(6) .+ x[14] - parameters.absorption_params.vmr["CO2"] = [a1; a2; a3] - model, lin = lin_model_from_parameters(parameters); + #dVMR_H2O = zeros(FT, 2, length(q)) # x[8], x[9] + #dVMR_CO2 = zeros(FT, 7, length(q)) # x[1], x[2], x[3], x[4], x[5], x[6], x[7] + + + parameters.absorption_params.vmr["H2O"] = [zeros(65).+ x[8]; + zeros(length(parameters.q)-65) .+ x[9]]; + #a1 = zeros(7) .+ x[12] + #a2 = zeros(7) .+ x[13] + #a3 = zeros(6) .+ x[14] + #H = 8.0 # km (atmospheric scale height) + #z = H*log.(psurf./parameters.p) + #prof = LogNormal(x[6], x[7]) + + parameters.absorption_params.vmr["CO2"] = (x[2].+zeros(length(parameters.q))) #+ + #(x[3] * exp.(-z./x[5])) + + #(x[4] * pdf.(prof, z)) #[a1; a2; a3] + + #dVMR_H2O[1,:] = [parameters.q[1:65] * 1.8; parameters.q[66:end] * 0.0;] # wrt x[7] + #dVMR_H2O[2,:] = [parameters.q[1:65] * 0.0; parameters.q[66:end] * 1.8;] # wrt x[8] + + #dVMR_CO2[1,:] = x[3] * exp.(-z./x[5]) * (-H/x[5]) * (1/x[1]); # wrt x[1] + #dVMR_CO2[2,:] = 1.0 .+ zeros(length(z)) # wrt x[2] + #dVMR_CO2[3,:] = exp.(-z./x[5]) # wrt x[3] + #dVMR_CO2[4,:] = pdf.(prof, z) # wrt x[4] + #dVMR_CO2[5,:] = x[3] * exp.(-z./x[5]) .* z./(x[5])^2 # wrt x[5] + #dVMR_CO2[6,:] = x[4] * pdf.(prof, z) .* (log.(z) .- x[6]) / x[7]^2 + #dVMR_CO2[7,:] = (x[4] * pdf.(prof, z) / x[7]) .* + # (((log.(z) - x[6]) / x[7]).^2 .- 1) + + model, lin = lin_model_from_parameters(parameters, x[1:9])#, dVMR_CO2, dVMR_H2O); RS_type = InelasticScattering.noRS( fscattRayl = [FT(1)], ϖ_Cabannes = [FT(1)], diff --git a/test/test_parameters/3BandParameters.yaml b/test/test_parameters/3BandParameters.yaml index eb0add75..58685a27 100644 --- a/test/test_parameters/3BandParameters.yaml +++ b/test/test_parameters/3BandParameters.yaml @@ -4,9 +4,9 @@ radiative_transfer: # Spectral bands (list of ν_start:ν_step:ν_end, in cm⁻¹) spec_bands: - - (1e7/773):0.02:(1e7/757) - - (1e7/1622):0.02:(1e7/1589) - - (1e7/2084):0.01:(1e7/2042) + - (1e7/773):0.05:(1e7/757) + - (1e7/1622):0.05:(1e7/1589) + - (1e7/2084):0.05:(1e7/2042) # Bidirectional Reflectance Distribution Function (BRDF) per band surface: - LambertianSurfaceScalar(0.2377) @@ -63,7 +63,7 @@ atmospheric_profile: 176.85, 236.64, 314.58, 418.87, 557.76, 735.00, 800.12, 849.00, 912.00, 980.00, 1005.0] # Reduce profile to n layers - profile_reduction: -1 + profile_reduction: 5 #-1 # ================================================================= # Absorption-Related Parameters (Optional) @@ -98,13 +98,20 @@ absorption: scattering: # List of scattering aerosols and their properties aerosols: - - τ_ref: 0.1 # Reference τ + - τ_ref: 0.1 #0.05 # Reference τ + μ: 0.1 #1.0 # effective radius (µm) + σ: 1.5 #1.2 # geometric standard deviation (µm) + nᵣ: 1.3 #1.6 # Real part of refractive index + nᵢ: 0.00000001 #0.0001 # Imag part of refractive index + z₀: 3.0 #2.0 #700.0 # Pressure peak (hPa) + σ₀: 0.6 #0.6 #50.0 # Pressure peak width (hPa) + - τ_ref: 0.05 # Reference τ μ: 1.0 # effective radius (µm) - σ: 1.5 # geometric standard deviation (µm) - nᵣ: 1.3 # Real part of refractive index + σ: 1.2 # geometric standard deviation (µm) + nᵣ: 1.6 # Real part of refractive index nᵢ: 0.00000001 # Imag part of refractive index - p₀: 700.0 # Pressure peak (hPa) - σp: 50.0 # Pressure peak width (hPa) + z₀: 2.0 #700.0 # Pressure peak (hPa) + σ₀: 0.6 #50.0 # Pressure peak width (hPa) # Maximum aerosol particle radius for quadrature points/weights (µm) r_max: 50.0 # Number of quadrature points for aerosol radius