Aerosol MMR set by condensation of requested species. Add coldtrap option to config. Error->warn in solver functions

Author: nichollsh

docs/src/reference/configuration.md

@@ -56,7 +56,7 @@ ratios (`vmr_dict` or `vmr_file`), or by metallicities (`metallicities`).
 | `metallicities`   | Dictionary of elemental **mass** abundance ratios relative to hydrogen. Must also set `p_surf`. |
 | `condensates   `  | List of volatiles which are allowed to condense. Can be used together with thermochemical equilibrium (`physics.chemistry = true`). |
 | `transparent   `  | Make the atmosphere transparent (see below). Replaces all of the above parameters in this table. |
-| `aerosols      `  | Dictionary of aerosols and their mass mixing ratios [kg/kg]. |
+| `aerosols      `  | Dictionary of aerosols, and their mass mixing ratios [kg/kg] or associated condensate species. |
 
 
 
@@ -104,6 +104,7 @@ Parameters that describe how the model should treat the physics.
 | `convection_crit` | Criterion for convective stability. Options: (s)chwarzschild, (l)edoux. |
 | `latent_heat   `  | Include vertical heat transport from condensation and evaporation (true/false). |
 | `rainout       `  | Enable condensation and evaporation of condensables aloft. Required for `latent_heat=true`. |
+| `coldtrap      `  | Enable cold-trapping effect on abundances aloft, so that VMR always decreases with height. |
 | `oceans        `  | Enable condensation and evaporation of condensables at the surface. |
 
 ## `[plots]`

res/config/default.toml

@@ -68,6 +68,7 @@ title = "Default"                       # Name for this configuration file
     conduction      = true                      # Include thermal conductive heat transport?
     oceans          = false                     # Model ocean formation - check saturation at surface
     rainout         = false                     # Model rainout - phase change impacts gas mixing ratios, not just energy fluxes
+    coldtrap        = true                      # Enable cold-trapping effect on abundances aloft
     latent_heat     = false                     # Include heat release from phase change
     convection_crit = "s"                       # Stability criterion for convection. Options: (s)chwarzschild or (l)edoux.
     convection      = true                      # Include heat transport by convection

res/config/physics/aerosols.toml

@@ -2,7 +2,7 @@
 title = "Aerosols radiative properties test"
 
 [planet]
-    tmp_surf        = 2000.0            # Surface temperature [kelvin].
+    tmp_surf        = 600.0            # Surface temperature [kelvin].
     instellation    = 44000.0           # Stellar flux at planet's orbital distance [W m-2].
     albedo_b        = 0.18              # Pseudo bond-albedo which downscales the stellar flux by 1-this_value.
     s0_fact         = 0.6652            # Stellar flux scale factor which accounts for planetary rotation (c.f. Cronin+13).
@@ -26,15 +26,15 @@ title = "Aerosols radiative properties test"
     io_dir          = "out/"                                            # Path for fast file operations (optional).
 
 [composition]
-    p_surf          = 270.0                     # Total surface pressure [bar].
+    p_surf          = 1000.0                     # Total surface pressure [bar].
     p_top           = 1e-5                      # Total top-of-atmosphere pressure [bar].
-    vmr_dict        = { H2O=0.05, CH4=0.1, CO=0.05, N2=0.8 }               # Volatile volume mixing ratios (=mole fractions).
-    condensates     = []                        # List of volatiles which are allowed to condense.
-    aerosols        = { "soot"=1e-4, "nitrate"=1e-2, "sulph"=1e-7 }            # List of aerosols and default MMRs
+    vmr_dict        = { H2=0.3, HCN=0.5, NH3=0.1, CO=0.05, H2O=0.05}               # Volatile volume mixing ratios (=mole fractions).
+    condensates     = ["HCN","NH3"]                        # List of volatiles which are allowed to condense.
+    aerosols        = { "soot"=1e-4, "nitrate"=1e-2, "biogenic"="HCN" }            # List of aerosols and default MMRs
 
 [execution]
     clean_output    = true                      # Clean the output folder at model startup.
-    verbosity       = 1                         # Log level (0: none, 1: normal, 2: debug)
+    verbosity       = 2                         # Log level (0: none, 1: normal, 2: debug)
     max_steps       = 20000                     # Maximum number of solver steps.
     max_runtime     = 400                       # Maximum wall-clock solver runtime [s].
     num_levels      = 200                       # Number of model levels.
@@ -44,7 +44,7 @@ title = "Aerosols radiative properties test"
     solution_type   = 0                         # Solution type (see wiki).
     solver          = ""                        # Solver to use (see wiki).
     dx_max          = 200.0                     # Maximum step size [Kelvin], when using nonlinear solvers
-    initial_state   = ["dry", "sat", "H2O"]     # Ordered list of requests describing the initial state of the atmosphere (see wiki).
+    initial_state   = ["dry"]     # Ordered list of requests describing the initial state of the atmosphere (see wiki).
     linesearch      = 0                         # Linesearch method to be used (0: None, 1: Golden, 2: Backtracking)
     easy_start      = false                     # Initially down-scale convective/condensation fluxes, if initial guess is poor.
     grey_start      = false                     # Initially solve for RCE using double-grey rad trans.
@@ -64,7 +64,8 @@ title = "Aerosols radiative properties test"
     sensible_heat   = false                     # Include sensible heat transport at the surface?
     conduction      = true                      # Include thermal conductive heat transport?
     oceans          = false                     # Model ocean formation - check saturation at surface
-    rainout         = false                     # Model rainout - phase change impacts gas mixing ratios, not just energy fluxes
+    rainout         = true                      # Model rainout - phase change impacts gas mixing ratios, not just energy fluxes
+    coldtrap        = true                      # Enable cold-trapping effect on abundances aloft
     latent_heat     = false                     # Include heat release from phase change
     convection_crit = "s"                       # Stability criterion for convection. Options: (s)chwarzschild or (l)edoux.
     convection      = true                      # Include heat transport by convection

src/AGNI.jl

@@ -267,6 +267,7 @@ module AGNI
         chem::Bool                          = false
         rainout::Bool                       = false
         oceans::Bool                        = false
+        coldtrap::Bool                      = true
         p_surf::Float64                     = 0.0
         p_top::Float64                      = 0.0
         pp_dict::Dict{String, Float64}      = Dict{String, Float64}()
@@ -297,6 +298,11 @@ module AGNI
             oceans  = Bool(cfg["physics"]["oceans"])
             condensates = cfg["composition"]["condensates"]
 
+            coldtrap = true
+            if haskey(cfg["physics"], "coldtrap")
+                coldtrap = Bool(cfg["physics"]["coldtrap"])
+            end
+
             comp_set_by::Int = 0
 
             # composition set by VMRs and Psurf
@@ -490,10 +496,10 @@ module AGNI
         end
 
         # Optional aerosol parametrization controls
-        aerosol_species::Dict{String, Float64} = Dict{String, Float64}()
+        aerosol_species::Dict = Dict()
         if haskey(cfg["composition"], "aerosols")
             for (k, v) in cfg["composition"]["aerosols"]
-                aerosol_species[string(k)] = Float64(v)
+                aerosol_species[string(k)] = v
             end
         end
 
@@ -518,6 +524,7 @@ module AGNI
 
                                 IO_DIR=io_dir,
                                 condensates=condensates,
+                                coldtrap=coldtrap,
                                 metallicities=metallicities,
                                 flag_gcontinuum   = cfg["physics"]["continua"],
                                 flag_rayleigh     = cfg["physics"]["rayleigh"],
@@ -641,6 +648,7 @@ module AGNI
                                 conv_rtol=conv_rtol,
                                 method=Int(method_idx),
                                 rainout=rainout,
+                                coldtrap=coldtrap,
                                 oceans=oceans,
                                 dx_max=Float64(cfg["execution"]["dx_max"]),
                                 ls_method=Int(cfg["execution"]["linesearch"]),

src/atmosphere.jl

@@ -290,6 +290,7 @@ module atmosphere
         aerosol_arr_l::Dict{String, Array{Float64,1}}  # Aerosol mass mixing ratio profiles [kg/kg]
         aerosol_arr_r::Dict{String, Array{Float64,1}}  # Aerosol particle size profiles [m]
         aerosol_val_r::Float64                         # Default particle size for aerosol species, if not specified in array
+        aerosol_setby::Dict{String, String}            # Dict of how each aerosol is set (e.g. "value", "S8", "H2O", etc.)
         aerosol_names::Array{String,1}               # Map SOCRATES index (int) to name (string)
         aerosol_relhumid::Float64                    # Mean relative humidity used by moist aerosol schemes [0,1]
         aerosol_phase_num::Int                       # Number of phase-function moments retained when averaging
@@ -440,7 +441,7 @@ module atmosphere
     - `flag_gcontinuum::Bool`           include generalised continuum absorption?
     - `flag_continuum::Bool`            include continuum absorption?
     - `flag_aerosol::Bool`              include aerosols?
-    - `aerosol_species::Dict`           aerosols MMR values to initialise profile with
+    - `aerosol_species::Dict`           aerosols MMR values or associated condensates
     - `flag_cloud::Bool`                include clouds?
     - `phs_timescale::Float64`          phase change timescale [s]
     - `evap_efficiency::Float64`        re-evaporatione efficiency compared to saturating amount
@@ -497,7 +498,7 @@ module atmosphere
                     flag_continuum::Bool =      false,
                     flag_aerosol::Bool =        false,
                     flag_cloud::Bool =          false,
-                    aerosol_species::Dict =     Dict{String, Float64}(),
+                    aerosol_species::Dict =     Dict{String, Union{Float64,String}}(),
 
                     phs_timescale::Float64 =    1e6,
                     evap_efficiency::Float64 =  0.05,
@@ -807,16 +808,33 @@ module atmosphere
         atmos.aerosol_val_r = 1.0e-5   # [INPUT] default particle size for aerosol species
         atmos.aerosol_arr_l = Dict{String, Array{Float64,1}}() # list of MMR profiles
         atmos.aerosol_arr_r = Dict{String, Array{Float64,1}}() # list of particle size profiles
+        atmos.aerosol_setby = Dict{String, String}() # dictionary of how each aerosol is set (e.g. "value", "S8", "H2O", etc.)
         atmos.aerosol_names = String[] # list of species names, in same order as spectral file
         for (k, v) in aerosol_species
+
+            # set to zero for now (true values will be set elsewhere)
             k = lowercase(k)
-            _check_range("Aerosol mass mixing ratio override for type $k", v; min=0.0) || return false
             atmos.aerosol_arr_l[k] = zeros(Float64, atmos.nlev_c)
             atmos.aerosol_arr_r[k] = zeros(Float64, atmos.nlev_c)
-            set_aerosol!(atmos,k , v)
+            push!(atmos.aerosol_names, UNSET_STR) # set in allocate!
 
-            # Store empty strings for now (set in allocate! function)
-            push!(atmos.aerosol_names, UNSET_STR)
+            # MMR profile tied to condensate or set by value?
+            if isa(v, String)
+                # interpret as condensate name
+                if !(v in condensates)
+                    @error "Aerosol '$k' is tied to '$v', but $v is not in the list of condensates"
+                    return false
+                end
+                @debug "Aerosol '$k' to be associated with species '$v'"
+                set_aerosol!(atmos, k, 0.0)
+                atmos.aerosol_setby[k] = v
+            else
+                # interpret as MMR value
+                v = Float64(v)
+                _check_range("Aerosol mass mixing ratio override for type $k", v; min=0.0) || return false
+                set_aerosol!(atmos, k, v)
+                atmos.aerosol_setby[k] = "value"
+            end
         end
 
         # Read VMRs
@@ -2819,6 +2837,25 @@ module atmosphere
         return aerosol_names
     end
 
+    """
+    **Calculate condensate mass mixing ratio at a specific layer.**
+
+    Arguments:
+    - `atmos::atmosphere.Atmos_t`   atmosphere struct instance to be used.
+    - `c::String`                   condensate species to calculate for (e.g. "H2O")
+    - `i::Int`                      layer index
+
+    Returns:
+    - `cond_mmr::Float64`           condensate `c` mass mixing ratio at layer `i`
+    """
+    function calc_cond_mmr(atmos::atmosphere.Atmos_t, c::String, i::Int64)::Float64
+        if atmos.cond_yield[c][i] > 0.0
+            return atmos.cond_yield[c][i]*atmos.cloud_alpha / atmos.layer_σ[i]
+        else
+            return 0.0
+        end
+    end
+
     """
     **Set water cloud profile based on saturation.**
 
@@ -2841,8 +2878,7 @@ module atmosphere
             # liquid water content (take ratio of mass surface densities [kg/m^2])
             if from_yield
                 # set by condensation yield
-                atmos.cloud_arr_l[i] = (atmos.cond_yield["H2O"][i]*atmos.cloud_alpha) /
-                                        atmos.layer_σ[i]
+                atmos.cloud_arr_l[i] = calc_cond_mmr(atmos, "H2O", i)
             else
                 # set by mask
                 atmos.cloud_arr_l[i] = atmos.gas_sat["H2O"][i] ? atmos.cloud_val_l : 0.0
@@ -2864,27 +2900,37 @@ module atmosphere
     end
 
     """
-    **Set aerosol profiles for a given species.**
+    **Set aerosol MMR profiles for a given species.**
 
     Arguments:
     - `atmos::atmosphere.Atmos_t`   the atmosphere struct instance to be used
     - `species::String`             the name of the aerosol species to set
-    - `mmr`                         the mixing ratio (profile or scalar) to set
+    - `mmr`                         the mixing ratio to assign (1D array, or Float, or species String)
 
     Optional arguments
     - `pmin::Float64`  the minimum pressure [Pa]
     - `pmax::Float64`  the maximum pressure [Pa]
 
     """
     function set_aerosol!(atmos::atmosphere.Atmos_t, species::String,
-                            mmr::Union{Array{Float64, 1}, Float64} = 0.0;
+                            mmr::Union{Array{Float64, 1}, Float64, String} = 0.0;
                             pmin::Float64 = 0.0, pmax::Float64 = 1e9)::Bool
 
+
+        # Reset
+        fill!(atmos.aerosol_arr_l[species], 0.0)
+        fill!(atmos.aerosol_arr_r[species], 0.0)
+
         # Mask for pressure range
         idx_mask = (atmos.p .>= pmin) .& (atmos.p .<= pmax)
 
         # Set mixing ratio profile for aerosol
-        if isa(mmr, Float64)
+        if isa(mmr, String)
+            # set by species mmr
+            for i in collect(1:atmos.nlev_c)[idx_mask]
+                atmos.aerosol_arr_l[species][i] = atmosphere.calc_cond_mmr(atmos, mmr, i)
+            end
+        elseif isa(mmr, Float64)
             # constant profile
             atmos.aerosol_arr_l[species][idx_mask] .= mmr
         else

src/chemistry.jl

@@ -256,7 +256,7 @@ module chemistry
 
                 # Apply cold trap, implicitly accounting for condensable not making it
                 #   this high up by vertical dynamical-transport processes.
-                if atmos.gas_vmr[c][i] > maxvmr[c]
+                if (atmos.gas_vmr[c][i] > maxvmr[c]) && atmos.coldtrap
                     atmos.gas_vmr[c][i] = maxvmr[c]
                     atmos.gas_sat[c][i] = true
                 end
@@ -388,6 +388,13 @@ module chemistry
             atmosphere.set_cloud!(atmos; from_yield=true)
         end
 
+        # Set aerosol profiles at levels where condensation occurs
+        for aer in atmos.aerosol_names
+            if atmos.aerosol_setby[aer] != "value"
+                atmosphere.set_aerosol!(atmos, aer, atmos.aerosol_setby[aer])
+            end
+        end
+
         return nothing
     end
 
@@ -706,6 +713,11 @@ module chemistry
         # recalculate layer properties
         atmosphere.calc_layer_props!(atmos)
 
+        # Warn on failure
+        if state > 0
+            @warn "FastChem internal failure; elements may not be conserved (state=$state)"
+        end
+
         # See docstring for return codes
         return state
     end
@@ -770,10 +782,10 @@ module chemistry
             end
         end
 
-        # keep cold trapping effect on abundances?
-        if !atmos.coldtrap
-            reset_to_chem!(atmos)
-        end
+        # keep cold trapping rainout effect on abundances?
+        # if !atmos.coldtrap
+        #     reset_to_chem!(atmos)
+        # end
 
         return state
     end

src/solver.jl

@@ -961,30 +961,30 @@ module solver
             rm(path_plt, force=true)
             rm(path_jac, force=true)
         elseif code == CODE_ITE
-            @error "    failure (maximum iterations)"
+            @warn "    failure (maximum iterations)"
             plot_step()
         elseif code == CODE_SIN
-            @error "    failure (singular jacobian)"
+            @warn "    failure (singular jacobian)"
             plotting.jacobian(b, path_jac)
             plot_step()
         elseif code == CODE_TIM
-            @error "    failure (maximum time)"
+            @warn "    failure (maximum time)"
         elseif code == CODE_NAN
-            @error "    failure (NaN values)"
+            @warn "    failure (NaN values)"
             plot_step()
         elseif code == CODE_CFG
-            @error "    failure (configuration)"
+            @warn "    failure (configuration)"
         elseif code == CODE_OBJ
-            @error "    failure (objective function)"
+            @warn "    failure (objective function)"
             plot_step()
         elseif code == CODE_STP
-            @error "    failure (other; last step not ok)"
+            @warn "    failure (other; last step not ok)"
             plot_step()
         elseif code == CODE_HYD
-            @error "    failure (hydrostatic integration)"
+            @warn "    failure (hydrostatic integration)"
             plot_step()
         else
-            @error "    failure (other)"
+            @warn "    failure (other)"
             plot_step()
         end