Add config option for rainout

CITATION.cff

@@ -5,7 +5,7 @@ authors:
   given-names: "Harrison"
   orcid: "https://orcid.org/0000-0002-8368-4641"
 title: "AGNI"
-version: 0.10.1
+version: 0.11.0
 doi: 10.xx/xx.xx
-date-released: 2024-09-20
+date-released: 2024-11-19
 url: "https://github.com/nichollsh/AGNI"

Project.toml

@@ -1,7 +1,7 @@
 name = "AGNI"
 uuid = "ede838c1-9ec3-4ebe-8ae8-da4091b3f21c"
 authors = ["Harrison Nicholls <harrison.nicholls@physics.ox.ac.uk>"]
-version = "0.10.1"
+version = "0.11.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

codemeta.json

@@ -12,12 +12,12 @@
   ],
   "identifier": "",
   "codeRepository": "https://github.com/nichollsh/AGNI",
-  "datePublished": "2024-11-14",
-  "dateModified": "2024-11-14",
+  "datePublished": "2024-11-19",
+  "dateModified": "2024-11-19",
   "dateCreated": "2024-09-13",
   "description": "A radiative-convective model for lava planet atmospheres.",
   "keywords": "physics, radiative transfer, exoplanets, astronomy, convection, radiation, planets, atmospheres",
   "license": "GPL v3.0",
   "title": "AGNI",
-  "version": "0.10.1"
+  "version": "0.11.0"
 }

res/config/55cnce_chem.toml

@@ -42,6 +42,7 @@ title = "Roughly 55 Cancri e @ fO2=IW"
     sensible_heat   = false
     latent_heat     = true
     convection      = true
+    rainout         = false
     solution_type   = 3
     solvers         = ["newton"]
     dx_max          = 400.0

res/config/condense.toml

@@ -19,7 +19,7 @@ title = "Condensation test"
     tmp_magma       = 1700.0
 
 [files]
-    input_sf        = "res/spectral_files/nogit/Dayspring48/Dayspring.sf"
+    input_sf        = "res/spectral_files/Dayspring/48/Dayspring.sf"
     input_star      = "res/stellar_spectra/sun.txt"
     output_dir      = "out/"
 
@@ -47,6 +47,7 @@ title = "Condensation test"
     sensible_heat   = true
     latent_heat     = true
     convection      = true
+    rainout         = true
     solution_type   = 3
     solvers         = ["newton"]
     dx_max          = 200.0

res/config/default.toml

@@ -51,6 +51,7 @@ title = "Default"                       # Name for this configuration file
     sensible_heat   = false                     # Include sensible heat transport at the surface?
     latent_heat     = false                     # Include heat release from phase change
     convection      = true                      # Include heat transport by convection
+    rainout         = false                     # Model rainout - phase change impacts gas mixing ratios, not just energy fluxes
     solution_type   = 0                         # Solution type (see wiki).
     solvers         = []                        # Ordered list of solvers to apply (see wiki).
     dx_max          = 200.0                     # Maximum step size [Kelvin], when using nonlinear solvers

res/config/hotdry.toml

@@ -46,6 +46,7 @@ title = "Hot and dry"
     sensible_heat   = true
     latent_heat     = false
     convection      = true
+    rainout         = false
     solution_type   = 3
     solvers         = ["newton"]
     dx_max          = 400.0

res/config/k2141b.toml

@@ -0,0 +1,68 @@
+# AGNI configuration file
+title = "K2-141 b (IW+4)"
+
+[planet]
+    tmp_surf        = 3000.0
+    instellation    = 3949800
+    albedo_b        = 0.0
+    s0_fact         = 0.25
+    zenith_angle    = 54.74
+    surface_material= "res/surface_albedos/lunar_marebasalt.dat"
+    albedo_s        = 0.3
+    radius          = 9.57e+6
+    gravity         = 22.125
+    flux_int        = 0.0
+    turb_coeff      = 0.001
+    wind_speed      = 2.0
+    skin_k          = 2.0
+    skin_d          = 0.01
+    tmp_magma       = 3000
+
+[files]
+    input_sf        = "/network/group/aopp/planetary/RTP035_NICHOLLS_PROTEUS/fwl_data/spectral_files/Honeyside/256/Honeyside.sf"
+    input_star      = "/network/group/aopp/planetary/RTP035_NICHOLLS_PROTEUS/fwl_data/stellar_spectra/Named/hd85512.txt"
+    output_dir      = "out/"
+
+[composition]
+    p_top           = 1e-5
+    p_dict          = {H2O = 1.248 , CO2 = 314.959 , N2  = 21.342 , H2  = 0.035 , CO  = 70.207 , CH4 = 0.000 , S2  = 0.958 , SO2 = 32.921}
+    include_all     = false
+    chemistry       = 0
+    condensates     = []
+
+[execution]
+    clean_output    = true
+    verbosity       = 1
+    max_steps       = 20000
+    max_runtime     = 400
+    num_levels      = 50
+    continua        = true
+    rayleigh        = true
+    cloud           = false
+    aerosol         = false
+    overlap_method  = "ee"
+    thermo_funct    = true
+    gravity_funct   = true
+    sensible_heat   = true
+    latent_heat     = true
+    convection      = true
+    rainout         = false
+    solution_type   = 3
+    solvers         = ["newton"]
+    dx_max          = 400.0
+    initial_state   = ["loglin", "1800"]
+    linesearch      = 1
+    easy_start      = false
+    converge_atol   = 1.0e-2
+    converge_rtol   = 1.0e-3
+
+[plots]
+    at_runtime      = true
+    temperature     = true
+    fluxes          = true
+    contribution    = true
+    emission        = true
+    albedo          = true
+    mixing_ratios   = true
+    animate         = true
+    height          = true

res/config/lbl_validation.toml

@@ -46,6 +46,7 @@ title = "Line-by-line validation tests"
     sensible_heat   = false
     latent_heat     = false
     convection      = false
+    rainout         = false
     solution_type   = 1
     solvers         = []
     dx_max          = 400.0

res/config/mnras_chem.toml

@@ -46,6 +46,7 @@ title = "TRAPPIST-1 c (IW+5) at PROTEUS solidification"
     sensible_heat   = true
     latent_heat     = true
     convection      = true
+    rainout         = false
     solution_type   = 2
     solvers         = ["newton"]
     dx_max          = 400.0

src/AGNI.jl

@@ -271,7 +271,7 @@ module AGNI
         end
         if chem_type in [1,2,3]
             if length(condensates)>0
-                @error "Misconfiguration: FastChem coupling incompatible with condensation"
+                @error "Misconfiguration: chemistry coupling is incompatible with condensation"
                 return false
             else
                 mkdir(dir_fastchem)
@@ -302,6 +302,7 @@ module AGNI
         conv_rtol::Float64     =        cfg["execution"]["converge_rtol"]
         max_steps::Int         =        cfg["execution"]["max_steps"]
         max_runtime::Float64   =        cfg["execution"]["max_runtime"]
+        rainout::Bool          =        cfg["execution"]["rainout"]
 
         #    plotting stuff
         plt_run::Bool          = cfg["plots"]["at_runtime"]
@@ -491,7 +492,7 @@ module AGNI
             if sol == "none"
                 energy.calc_fluxes!(atmos, incl_latent,
                                     incl_convect, incl_sens, incl_conduct,
-                                    calc_cf=plt_cff)
+                                    calc_cf=plt_cff, rainout=rainout)
                 @info "    done"
 
             # Nonlinear solver
@@ -506,7 +507,7 @@ module AGNI
                                     sens_heat=incl_sens, max_steps=max_steps,
                                     max_runtime=max_runtime,
                                     conv_atol=conv_atol, conv_rtol=conv_rtol,
-                                    method=method_idx,
+                                    method=method_idx, rainout=rainout,
                                     dx_max=dx_max, ls_method=linesearch,
                                     easy_start=easy_start,
                                     modplot=modplot,save_frames=plt_ani)

src/atmosphere.jl

@@ -326,7 +326,7 @@ module atmosphere
         @info  "Setting-up a new atmosphere struct"
 
         # Code versions
-        atmos.AGNI_VERSION = "0.10.1"
+        atmos.AGNI_VERSION = "0.11.0"
         atmos.SOCRATES_VERSION = readchomp(joinpath(ENV["RAD_DIR"],"version"))
         @debug "AGNI VERSION = "*atmos.AGNI_VERSION
         @debug "Using SOCRATES at $(ENV["RAD_DIR"])"

src/energy.jl

@@ -333,7 +333,7 @@ module energy
     - `pmin::Float64`           pressure [bar] below which convection is disabled
     - `mltype::Int`             mixing length value (1: scale height, 2: asymptotic)
     """
-    function mlt!(atmos::atmosphere.Atmos_t; pmin::Float64=1.0e-4, mltype::Int=2)
+    function convection!(atmos::atmosphere.Atmos_t; pmin::Float64=1.0e-4, mltype::Int=2)
 
         pmin *= 1.0e5 # convert bar to Pa
 
@@ -602,7 +602,7 @@ module energy
             @turbo @. atmos.flux_l += atmos.phs_wrk_fl
 
             # calculate mask
-            @. atmos.mask_l = (abs(atmos.flux_l) > 1.0e-10)
+            @. atmos.mask_l = (abs(atmos.flux_l) > 1.0e-30)
 
         end # go to next condensable
 
@@ -652,6 +652,15 @@ module energy
         fill!(atmos.flux_tot, 0.0)
     end
 
+    """
+    **Reset mixing ratios to their original values**
+    """
+    function restore_composition!(atmos::atmosphere.Atmos_t)
+        for g in atmos.gas_names
+            @turbo @. atmos.gas_vmr[g] = atmos.gas_ovmr[g]
+        end
+    end
+
     """
     **Calculate total flux at each level.**
 
@@ -664,52 +673,52 @@ module energy
     - `convect_sf::Float64`             scale factor applied to convection fluxes
     - `latent_sf::Float64`              scale factor applied to phase change fluxes
     - `calc_cf::Bool`                   calculate LW contribution function?
-    - `reset_vmrs::Bool`                reset VMRs to dry values before radtrans and MLT
+    - `rainout::Bool`                   allow rainout ( do not reset VMRs to dry values )
     """
     function calc_fluxes!(atmos::atmosphere.Atmos_t,
                           latent::Bool, convect::Bool, sens_heat::Bool, conduct::Bool;
                           convect_sf::Float64=1.0, latent_sf::Float64=1.0,
-                          calc_cf::Bool=false, reset_vmrs::Bool=true)
+                          calc_cf::Bool=false, rainout::Bool=true)
 
         # Reset fluxes
         reset_fluxes!(atmos)
 
         # +Condensation and evaporation
         if atmos.condense_any && latent
 
+            # Restore mixing ratios
+            restore_composition!(atmos)
             atmosphere.calc_layer_props!(atmos)
 
             # Handle rainout
             atmosphere.handle_saturation!(atmos)
 
             # Calculate latent heat flux
-            energy.condense_diffuse!(atmos)
+            condense_diffuse!(atmos)
 
             # Modulate?
             atmos.flux_l *= latent_sf
 
             # Add to total flux
             @turbo @. atmos.flux_tot += atmos.flux_l
 
-            # Restore mixing ratios
-            if reset_vmrs
-                for g in atmos.gas_names
-                    @turbo @. atmos.gas_vmr[g] = atmos.gas_ovmr[g]
-                end
+            # Restore mixing ratios - do not allow rainout
+            if !rainout
+                restore_composition!(atmos)
             end
         end
         # Calculate layer properties
         atmosphere.calc_layer_props!(atmos)
 
         # +Radiation
-        energy.radtrans!(atmos, true, calc_cf=calc_cf)
-        energy.radtrans!(atmos, false)
+        radtrans!(atmos, true, calc_cf=calc_cf)
+        radtrans!(atmos, false)
         @turbo @. atmos.flux_tot += atmos.flux_n
 
         # +Dry convection
         if convect
             # Calc flux
-            energy.mlt!(atmos)
+            convection!(atmos)
 
             # Modulate?
             atmos.flux_cdry *= convect_sf
@@ -720,13 +729,13 @@ module energy
 
         # +Surface turbulence
         if sens_heat
-            energy.sensible!(atmos)
+            sensible!(atmos)
             atmos.flux_tot[end] += atmos.flux_sens
         end
 
         # +Conduction
         if conduct
-            energy.conduct!(atmos)
+            conduct!(atmos)
             @turbo @. atmos.flux_tot += atmos.flux_cdct
         end
 

src/solver.jl

@@ -77,7 +77,7 @@ module solver
     - `sens_heat::Bool`                 include sensible heating at the surface
     - `conduct::Bool`                   include conductive heat transport within the atmosphere
     - `latent::Bool`                    include latent heat exchange (condensation/evaporation)
-    - `reset_vmrs::Bool`                reset atmosphere to a well-mixed composition after the phase change calculation
+    - `rainout::Bool`                   allow rainout (phase change impacts mixing ratios, not just energy fluxes)
     - `dx_max::Float64`                 maximum step size [K]
     - `max_steps::Int`                  maximum number of solver steps
     - `max_runtime::Float64`            maximum runtime in wall-clock seconds
@@ -104,8 +104,7 @@ module solver
                             sol_type::Int=1,
                             chem_type::Int=0,
                             convect::Bool=true, sens_heat::Bool=true,
-                            conduct::Bool=true, latent::Bool=true,
-                            reset_vmrs::Bool=true,
+                            conduct::Bool=true, latent::Bool=true, rainout::Bool=true,
                             dx_max::Float64=400.0,
                             max_steps::Int=400, max_runtime::Float64=900.0,
                             fdw::Float64=3.0e-5, fdc::Bool=true, fdo::Int=2,
@@ -244,7 +243,7 @@ module solver
             energy.calc_fluxes!(atmos,
                                 latent, convect, sens_heat, conduct,
                                 convect_sf=easy_sf, latent_sf=easy_sf,
-                                reset_vmrs=reset_vmrs)
+                                rainout=rainout)
 
             # Energy divergence term
             @turbo @. atmos.flux_dif -= atmos.ediv_add