Skip to content

Commit

Permalink
fixing tests + more cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@trontrytel
trontrytel committed Aug 13, 2023
1 parent eb68a75 commit 6dcb4c2
Show file tree
Hide file tree
Showing 5 changed files with 88 additions and 97 deletions.
72 changes: 69 additions & 3 deletions src/Common.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,13 @@ import Thermodynamics as TD
import ..Parameters as CMP
const APS = CMP.AbstractCloudMicrophysicsParameters

import ..CommonTypes as CT

export G_func
export H2SO4_soln_saturation_vapor_pressure
export ABIFM_Delta_a_w
export Chen2022_vel_add
export Chen2022_vel_coeffs

"""
G_func(param_set, T, Liquid())
Expand Down Expand Up @@ -167,6 +170,70 @@ function Delta_a_w(prs::APS, x::FT, T::FT) where {FT <: Real}
return min(Δa_w, FT(1))
end

"""
Chen2022_vel_coeffs(prs, precip_type, ρ)
- prs - set with free parameters
- precip_type - type for ice, rain or snow
- ρ - air density
Returns the coefficients from Appendix B in Chen et al 2022
DOI: 10.1016/j.atmosres.2022.106171
"""
function Chen2022_vel_coeffs(prs::APS, ::CT.RainType, ρ::FT) where {FT <: Real}

ρ0::FT = CMP.q_coeff_rain_Ch2022(prs)
a1::FT = CMP.a1_coeff_rain_Ch2022(prs)
a2::FT = CMP.a2_coeff_rain_Ch2022(prs)
a3::FT = CMP.a3_coeff_rain_Ch2022(prs)
a3_pow::FT = CMP.a3_pow_coeff_rain_Ch2022(prs)
b1::FT = CMP.b1_coeff_rain_Ch2022(prs)
b2::FT = CMP.b2_coeff_rain_Ch2022(prs)
b3::FT = CMP.b3_coeff_rain_Ch2022(prs)
b_ρ::FT = CMP.b_rho_coeff_rain_Ch2022(prs)
c1::FT = CMP.c1_coeff_rain_Ch2022(prs)
c2::FT = CMP.c2_coeff_rain_Ch2022(prs)
c3::FT = CMP.c3_coeff_rain_Ch2022(prs)

q = exp(ρ0 * ρ)
ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow)
bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ)
ci = (c1, c2, c3)

# unit conversions
aiu = ai .* 1000 .^ bi
ciu = ci .* 1000

return(aiu, bi, ciu)
end
function Chen2022_vel_coeffs(
prs::APS,
::Union{CT.IceType, CT.SnowType},
ρ::FT,
) where {FT <: Real}

ρ_i::FT = CMP.ρ_cloud_ice(prs)

# TODO - move to CLIMAParameters
As = FT(0.00174079) * (log(ρ_i))^2 − FT(0.0378769) * log(ρ_i) - FT(0.263503)
Bs = FT(1) / (FT(0.575231) + FT(0.0909307) * log(ρ_i) + FT(0.515579) / sqrt(ρ_i))
Cs = FT(-0.345387) + FT(0.177362) * exp(FT(-0.000427794) * ρ_i) +
FT(0.00419647) * sqrt(ρ_i)
Es = FT(-0.156593) - FT(0.0189334) * (log(ρ_i))^2 + FT(0.1377817) * sqrt(ρ_i)
Fs = -exp(FT(-3.35641) - FT(0.0156199) * (log(ρ_i))^2 + FT(0.765337) * log(ρ_i))
Gs = FT(1) /
(FT(-0.0309715) + FT(1.55054) / (log(ρ_i)) - FT(0.518349) * log(ρ_i) / ρ_i)

ai = (Es * ρ^As, Fs * ρ^As)
bi = (Bs + ρ * Cs, Bs + ρ * Cs)
ci = (FT(0), Gs)
# unit conversions
aiu = ai .* 1000 .^ bi
ciu = ci .* 1000

return(aiu, bi, ciu)
end

"""
Chen2022_vel_add(a, b, c, λ, k)
Expand All @@ -180,8 +247,7 @@ We are assuming exponential size distribution and hence μ=0.
"""
function Chen2022_vel_add(a::FT, b::FT, c::FT, λ::FT, k::Int) where {FT <: Real}
μ = 0 # Exponantial instaed of gamma distribution
δ = FT(μ + k + 1)
return a * λ^δ * SF.gamma(b + δ) / (λ + c)^(b .+ δ) / SF.gamma(δ)
δ = μ + k + 1
return a * λ^δ * SF.gamma(b + δ) / (λ + c)^(b + δ) / SF.gamma(δ)
end

end
67 changes: 7 additions & 60 deletions src/Microphysics1M.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,28 +219,9 @@ function terminal_velocity(
) where {FT <: Real}
fall_w = FT(0)
if q_ > FT(0)
# TODO - AJ - try to unify rain and ice
# coefficients from Table B1 from Chen et. al. 2022
ρ0::FT = CMP.q_coeff_rain_Ch2022(prs)
a1::FT = CMP.a1_coeff_rain_Ch2022(prs)
a2::FT = CMP.a2_coeff_rain_Ch2022(prs)
a3::FT = CMP.a3_coeff_rain_Ch2022(prs)
a3_pow::FT = CMP.a3_pow_coeff_rain_Ch2022(prs)
b1::FT = CMP.b1_coeff_rain_Ch2022(prs)
b2::FT = CMP.b2_coeff_rain_Ch2022(prs)
b3::FT = CMP.b3_coeff_rain_Ch2022(prs)
b_ρ::FT = CMP.b_rho_coeff_rain_Ch2022(prs)
c1::FT = CMP.c1_coeff_rain_Ch2022(prs)
c2::FT = CMP.c2_coeff_rain_Ch2022(prs)
c3::FT = CMP.c3_coeff_rain_Ch2022(prs)
q = exp(ρ0 * ρ)
ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow)
bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ)
ci = (c1, c2, c3)
#unit conversions
aiu = ai .* 1000 .^ bi
ciu = ci .* 1000

# coefficients from Table B1 from Chen et. al. 2022
aiu, bi, ciu = CO.Chen2022_vel_coeffs(prs, precip, ρ)
# size distribution parameter
_λ::FT = lambda(prs, precip, q_, ρ)

Expand All @@ -262,29 +243,11 @@ function terminal_velocity(
if q_ > FT(0)

ρ_i::FT = CMP.ρ_cloud_ice(prs)

# TODO - move to CLIMAParameters
As = FT(0.00174079 * (log(ρ_i))^2 − 0.0378769 * log(ρ_i) - 0.263503)
Bs = FT((0.575231 + 0.0909307 * log(ρ_i) + 0.515579 / sqrt(ρ_i))^(-1))
Cs = FT(
-0.345387 +
0.177362 * exp(-0.000427794 * ρ_i) +
0.00419647 * sqrt(ρ_i))
Es = FT(-0.156593 - 0.0189334 * (log(ρ_i))^2 + 0.1377817 * sqrt(ρ_i))
Fs = FT(-exp(-3.35641 - 0.0156199 * (log(ρ_i))^2 + 0.765337 * log(ρ_i)))
Gs = FT(
(-0.0309715 + 1.55054 / (log(ρ_i)) - 0.518349 * log(ρ_i) / ρ_i)^(-1))

ai = (Es * ρ^As, Fs * ρ^As)
bi = (Bs + ρ * Cs, Bs + ρ * Cs)
ci = (FT(0), Gs)
#unit conversions
aiu = ai .* 1000 .^ bi
ciu = ci .* 1000

# size distribution parameter
_λ::FT = lambda(prs, precip, q_, ρ)

# coefficients from Appendix B from Chen et. al. 2022
aiu, bi, ciu = CO.Chen2022_vel_coeffs(prs, precip, ρ)

# eq 20 from Chen et al 2022
fall_w = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, _λ, 3))
# It should be ϕ^κ * fall_w, but for rain drops ϕ = 1 and κ = 0
Expand Down Expand Up @@ -312,24 +275,8 @@ function terminal_velocity(

ρ_i::FT = CMP.ρ_cloud_ice(prs)

# TODO - move to CLIMAParameters
As = 0.00174079 * (log(ρ_i))^2 − 0.0378769 * log(ρ_i) - 0.263503
Bs = (0.575231 + 0.0909307 * log(ρ_i) + 0.515579 / sqrt(ρ_i))^(-1)
Cs =
-0.345387 +
0.177362 * exp(-0.000427794 * ρ_i) +
0.00419647 * sqrt(ρ_i)
Es = -0.156593 - 0.0189334 * (log(ρ_i))^2 + 0.1377817 * sqrt(ρ_i)
Fs = -exp(-3.35641 - 0.0156199 * (log(ρ_i))^2 + 0.765337 * log(ρ_i))
Gs =
(-0.0309715 + 1.55054 / (log(ρ_i)) - 0.518349 * log(ρ_i) / ρ_i)^(-1)

ai = (Es * ρ^As, Fs * ρ^As)
bi = (Bs + ρ * Cs, Bs + ρ * Cs)
ci = (0, Gs)
#unit conversions
aiu = ai .* 1000 .^ bi
ciu = ci .* 1000
# coefficients from Appendix B from Chen et. al. 2022
aiu, bi, ciu = CO.Chen2022_vel_coeffs(prs, precip, ρ)

κ = FT(-1 / 3) #oblate
k = 3 # mass weighted
Expand Down
23 changes: 1 addition & 22 deletions src/Microphysics2M.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,7 +398,6 @@ function rain_terminal_velocity(
ρ::FT,
N_rai::FT,
) where {FT <: Real}

if q_rai < eps(FT)
return (FT(0), FT(0))
end
Expand All @@ -425,28 +424,8 @@ function rain_terminal_velocity(
if q_rai < eps(FT)
return (FT(0), FT(0))
end

# coefficients from Table B1 from Chen et. al. 2022
ρ0::FT = CMP.q_coeff_rain_Ch2022(param_set)
a1::FT = CMP.a1_coeff_rain_Ch2022(param_set)
a2::FT = CMP.a2_coeff_rain_Ch2022(param_set)
a3::FT = CMP.a3_coeff_rain_Ch2022(param_set)
a3_pow::FT = CMP.a3_pow_coeff_rain_Ch2022(param_set)
b1::FT = CMP.b1_coeff_rain_Ch2022(param_set)
b2::FT = CMP.b2_coeff_rain_Ch2022(param_set)
b3::FT = CMP.b3_coeff_rain_Ch2022(param_set)
b_ρ::FT = CMP.b_rho_coeff_rain_Ch2022(param_set)
c1::FT = CMP.c1_coeff_rain_Ch2022(param_set)
c2::FT = CMP.c2_coeff_rain_Ch2022(param_set)
c3::FT = CMP.c3_coeff_rain_Ch2022(param_set)
q = exp(ρ0 * ρ)
ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow)
bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ)
ci = (c1, c2, c3)
#unit conversions
aiu = ai .* 1000 .^ bi
ciu = ci .* 1000

aiu, bi, ciu = CO.Chen2022_vel_coeffs(param_set, CT.RainType(), ρ)
# size distribution parameter
λ = raindrops_limited_vars(param_set, q_rai, ρ, N_rai).λr

Expand Down
3 changes: 1 addition & 2 deletions test/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,5 @@ CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"

[compat]
CLIMAParameters = ">=0.7.5"
CLIMAParameters = ">=0.7.8"
20 changes: 10 additions & 10 deletions test/microphysics_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ const B1994 = CMT.B1994Type()
const TC1980 = CMT.TC1980Type()
const LD2004 = CMT.LD2004Type()
const Blk1MVel = CMT.Blk1MVelType()
const SB2006Vel = CMT.Chen2022Type()
const SB2006Vel = CMT.SB2006VelType()
const Ch2022 = CMT.Chen2022Type()

@info "Microphysics Tests"
Expand Down Expand Up @@ -117,7 +117,7 @@ function test_microphysics(FT)
v_bigger = CM1.terminal_velocity(prs, rain, Ch2022, ρ, q_rai * 2)

#test
TT.@test vt_rai ≈ 3.0721397260869705 rtol = 1e-6
TT.@test vt_rai ≈ 3.0721397260869705 rtol = 2e-6
TT.@test CM1.terminal_velocity(prs, rain, Ch2022, ρ, FT(0))[1] ≈ 0 atol =
eps(FT)
TT.@test v_bigger > vt_rai
Expand All @@ -133,12 +133,13 @@ function test_microphysics(FT)
vt_ice = CM1.terminal_velocity(prs, ice, Ch2022, ρ, q_ice)
v_bigger = CM1.terminal_velocity(prs, ice, Ch2022, ρ, q_ice * 2)

# TODO -TODO!!!
# TODO-term_vel - where does the 4.045213370181404 come from?
#test
println("vt_ice = ", vt_ice)
#TT.@test vt_ice ≈ 4.045213370181404 rtol = 1e-6
#TT.@test CM1.terminal_velocity(prs, ice, Ch2022, ρ, FT(0)) ≈ 0 atol =
# eps(FT)
#TT.@test v_bigger > vt_ice
TT.@test CM1.terminal_velocity(prs, ice, Ch2022, ρ, FT(0)) ≈ 0 atol =
eps(FT)
TT.@test v_bigger > vt_ice
end

TT.@testset "1M_microphysics - Chen 2022 snow terminal velocity" begin
Expand All @@ -152,7 +153,7 @@ function test_microphysics(FT)
v_bigger = CM1.terminal_velocity(prs, snow, Ch2022, ρ, q_sno * 2)

#test
TT.@test vt_sno ≈ 1.4722373984934243 rtol = 1e-6
TT.@test vt_sno ≈ 1.4722373984934243 rtol = 2e-6
TT.@test CM1.terminal_velocity(prs, snow, Ch2022, ρ, FT(0)) ≈ 0 atol =
eps(FT)
TT.@test v_bigger > vt_sno
Expand Down Expand Up @@ -865,9 +866,8 @@ function test_microphysics(FT)

#test
TT.@test vt_rai isa Tuple
# TODO!!!
#TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6
#TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6
TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6
TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6
TT.@test CM2.rain_terminal_velocity(
prs,
SB2006,
Expand Down

0 comments on commit 6dcb4c2

Please sign in to comment.