Version 0.8.4. New plot for pressure vs height (new cfg option). Pres…

…sure grid uses a smaller topmost level to avoid issues with extrapolation with low opacity compositions. New tests for making plots. Flag for updating (or not) the mixing ratios for radtrans post-saturation check. Fixed loop instability from SIMD when reading SOCRATES radout. Ensures that input stellar spectrum is strictly ascending in wavelength. Minor changes to notebooks

CITATION.cff

@@ -5,7 +5,7 @@ authors:
   given-names: "Harrison"
   orcid: "https://orcid.org/0000-0002-8368-4641"
 title: "AGNI"
-version: 0.8.3
+version: 0.8.4
 doi: 10.xx/xx.xx
 date-released: 2024-09-20
 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.8.3"
+version = "0.8.4"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

codemeta.json

@@ -19,5 +19,5 @@
   "keywords": "physics, radiative transfer, exoplanets, astronomy, convection, radiation, planets, atmospheres",
   "license": "GPL v3.0",
   "title": "AGNI",
-  "version": "0.8.3"
+  "version": "0.8.4"
 }

res/config/55cnce_chem.toml

@@ -59,4 +59,4 @@ title = "Roughly 55 Cancri e @ fO2=IW"
     albedo          = true
     mixing_ratios   = true
     animate         = true
-
+    height          = true

res/config/L-98-59d.toml

@@ -60,4 +60,4 @@ title = "L 98-59 d"
     albedo          = true
     mixing_ratios   = true
     animate         = true
-
+    height          = true

res/config/condense.toml

@@ -65,4 +65,4 @@ title = "Condensation test"
     albedo          = true
     mixing_ratios   = true
     animate         = true
-
+    height          = true

res/config/default.toml

@@ -68,4 +68,4 @@ title = "Default"                       # Name for this configuration file
     albedo          = true                      # Plot spectral albedo?
     mixing_ratios   = true                      # Plot mixing ratios?
     animate         = false                     # Make animation from runtime plots?
-
+    height          = true                      # Plot hydrostatic solution for height?

res/config/hotdry.toml

@@ -64,4 +64,4 @@ title = "Hot and dry"
     albedo          = true
     mixing_ratios   = true
     animate         = true
-
+    height          = true

src/AGNI.jl

@@ -282,35 +282,36 @@ module AGNI
         spfile_name::String    =        cfg["files" ]["input_sf"]
         star_file::String      =        cfg["files" ]["input_star"]
         nlev_centre::Int       =        cfg["execution"]["num_levels"]
-        flag_cnt::Bool         =        cfg["execution" ]["continua"]
-        flag_ray::Bool         =        cfg["execution" ]["rayleigh"]
-        flag_cld::Bool         =        cfg["execution" ]["cloud"]
-        flag_aer::Bool         =        cfg["execution" ]["aerosol"]
-        overlap::Int           =        cfg["execution" ]["overlap_method"]
-        thermo_funct::Bool     =        cfg["execution" ]["thermo_funct"]
-        incl_convect::Bool     =        cfg["execution" ]["convection"]
-        incl_sens::Bool        =        cfg["execution" ]["sensible_heat"]
-        incl_latent::Bool      =        cfg["execution" ]["latent_heat"]
-        sol_type::Int          =        cfg["execution" ]["solution_type"]
-        solvers_cmd::Array{String,1} =  cfg["execution" ]["solvers"]
-        initial_req::Array{String,1} =  cfg["execution" ]["initial_state"]
-        dx_max::Float64        =        cfg["execution" ]["dx_max"]
-        linesearch::Int        =        cfg["execution" ]["linesearch"]
-        easy_start::Bool       =        cfg["execution" ]["easy_start"]
-        conv_atol::Float64     =        cfg["execution" ]["converge_atol"]
-        conv_rtol::Float64     =        cfg["execution" ]["converge_rtol"]
-        max_steps::Int         =        cfg["execution" ]["max_steps"]
-        max_runtime::Float64   =        cfg["execution" ]["max_runtime"]
+        flag_cnt::Bool         =        cfg["execution"]["continua"]
+        flag_ray::Bool         =        cfg["execution"]["rayleigh"]
+        flag_cld::Bool         =        cfg["execution"]["cloud"]
+        flag_aer::Bool         =        cfg["execution"]["aerosol"]
+        overlap::Int           =        cfg["execution"]["overlap_method"]
+        thermo_funct::Bool     =        cfg["execution"]["thermo_funct"]
+        incl_convect::Bool     =        cfg["execution"]["convection"]
+        incl_sens::Bool        =        cfg["execution"]["sensible_heat"]
+        incl_latent::Bool      =        cfg["execution"]["latent_heat"]
+        sol_type::Int          =        cfg["execution"]["solution_type"]
+        solvers_cmd::Array{String,1} =  cfg["execution"]["solvers"]
+        initial_req::Array{String,1} =  cfg["execution"]["initial_state"]
+        dx_max::Float64        =        cfg["execution"]["dx_max"]
+        linesearch::Int        =        cfg["execution"]["linesearch"]
+        easy_start::Bool       =        cfg["execution"]["easy_start"]
+        conv_atol::Float64     =        cfg["execution"]["converge_atol"]
+        conv_rtol::Float64     =        cfg["execution"]["converge_rtol"]
+        max_steps::Int         =        cfg["execution"]["max_steps"]
+        max_runtime::Float64   =        cfg["execution"]["max_runtime"]
 
         #    plotting stuff
-        plt_run::Bool          = cfg["plots"     ]["at_runtime"]
-        plt_tmp::Bool          = cfg["plots"     ]["temperature"]
-        plt_flx::Bool          = cfg["plots"     ]["fluxes"]
-        plt_cff::Bool          = cfg["plots"     ]["contribution"]
-        plt_ems::Bool          = cfg["plots"     ]["emission"]
-        plt_alb::Bool          = cfg["plots"     ]["albedo"]
-        plt_vmr::Bool          = cfg["plots"     ]["mixing_ratios"]
-        plt_ani::Bool          = cfg["plots"     ]["animate"]
+        plt_run::Bool          = cfg["plots"]["at_runtime"]
+        plt_tmp::Bool          = cfg["plots"]["temperature"]
+        plt_flx::Bool          = cfg["plots"]["fluxes"]
+        plt_cff::Bool          = cfg["plots"]["contribution"]
+        plt_ems::Bool          = cfg["plots"]["emission"]
+        plt_alb::Bool          = cfg["plots"]["albedo"]
+        plt_vmr::Bool          = cfg["plots"]["mixing_ratios"]
+        plt_hei::Bool          = cfg["plots"]["height"]
+        plt_ani::Bool          = cfg["plots"]["animate"]
         plt_ani = plt_ani && plt_tmp
 
         # Read OPTIONAL configuration options from dict
@@ -537,6 +538,7 @@ module AGNI
         plt_flx && plotting.plot_fluxes(atmos,   joinpath(atmos.OUT_DIR,"plot_fluxes.png"), incl_mlt=incl_convect, incl_eff=(sol_type==3), incl_cdct=incl_conduct, incl_latent=incl_latent)
         plt_ems && plotting.plot_emission(atmos, joinpath(atmos.OUT_DIR,"plot_emission.png"))
         plt_alb && plotting.plot_albedo(atmos,   joinpath(atmos.OUT_DIR,"plot_albedo.png"))
+        plt_hei && plotting.plot_height(atmos,   joinpath(atmos.OUT_DIR,"plot_height.png"))
 
         # Deallocate atmosphere
         @info "Deallocating memory"

src/atmosphere.jl

@@ -318,7 +318,7 @@ module atmosphere
         end
 
         # Code versions
-        atmos.AGNI_VERSION = "0.8.3"
+        atmos.AGNI_VERSION = "0.8.4"
         atmos.SOCRATES_VERSION = readchomp(joinpath(ENV["RAD_DIR"],"version"))
         @debug "AGNI VERSION = $(atmos.AGNI_VERSION)"
         @debug "Using SOCRATES at $(ENV["RAD_DIR"])"
@@ -423,8 +423,8 @@ module atmosphere
         atmos.cloud_arr_f   = zeros(Float64, atmos.nlev_c)
 
         # Phase change timescales [seconds]
-        atmos.phs_tau_mix = 1.0e5   # mixed composition case
-        atmos.phs_tau_sgl = 1.0e5   # single gas case
+        atmos.phs_tau_mix = 3.0e4   # mixed composition case
+        atmos.phs_tau_sgl = 3.0e4   # single gas case
 
         # Hardcoded cloud properties
         atmos.cond_alpha    = 0.0     # 0% of condensate is retained (i.e. complete rainout)
@@ -612,7 +612,7 @@ module atmosphere
 
         # store condensates
         for c in condensates
-            if !atmos.gas_dat[c].stub && !atmos.gas_dat[c].no_sat
+            if !atmos.gas_dat[c].stub && !atmos.gas_dat[c].no_sat && !(c == "H2")
                 push!(atmos.condensates, c)
             end
         end
@@ -850,8 +850,6 @@ module atmosphere
         atmos.p  = zeros(Float64, atmos.nlev_c)
         atmos.pl = zeros(Float64, atmos.nlev_l)
 
-        # First, assign log10'd values...
-
         # Top and bottom boundaries
         atmos.pl[end] = log10(atmos.p_boa)
         atmos.pl[1]   = log10(atmos.p_toa)
@@ -874,7 +872,11 @@ module atmosphere
         # Set cell-centres at midpoint of cell-edges
         atmos.p[1:end] .= 0.5 .* (atmos.pl[1:end-1] .+ atmos.pl[2:end])
 
-        # Finally, convert arrays to 'real' pressure units
+        # Shrink top-most layer to avoid doing too much extrapolation
+        p_fact = 0.8
+        atmos.p[1] = atmos.pl[1]*p_fact + atmos.p[1]*(1-p_fact)
+
+        # Finally, convert arrays to actual pressure units [Pa]
         @turbo @. atmos.p  = 10.0 ^ atmos.p
         @turbo @. atmos.pl = 10.0 ^ atmos.pl
 

src/energy.jl

@@ -227,7 +227,7 @@ module energy
         if lw
             # LW case
             for lv in 1:atmos.nlev_l      # sum over levels
-                @inbounds for ba in 1:atmos.dimen.nd_channel  # sum over bands
+                for ba in 1:atmos.dimen.nd_channel  # sum over bands
                     idx = lv+(ba-1)*atmos.nlev_l
                     atmos.band_d_lw[lv,ba] = max(0.0, atmos.radout.flux_down[idx])
                     atmos.band_u_lw[lv,ba] = max(0.0, atmos.radout.flux_up[idx])
@@ -242,7 +242,7 @@ module energy
             fill!(atmos.contfunc_norm,0.0)
             if calc_cf
                 cf_max::Float64 = maximum(atmos.radout.contrib_funcf_band[1,:,:])
-                @inbounds for ba in 1:atmos.dimen.nd_channel
+                for ba in 1:atmos.dimen.nd_channel
                     for lv in 1:atmos.nlev_c
                         atmos.contfunc_norm[lv,ba] =
                                             atmos.radout.contrib_funcf_band[1,lv,ba]/cf_max
@@ -254,7 +254,7 @@ module energy
         else
             # SW case
             for lv in 1:atmos.nlev_l                # sum over levels
-                @inbounds for ba in 1:atmos.dimen.nd_channel  # sum over bands
+                for ba in 1:atmos.dimen.nd_channel  # sum over bands
                     idx = lv+(ba-1)*atmos.nlev_l
                     atmos.band_d_sw[lv,ba] = max(0.0,atmos.radout.flux_down[idx])
                     atmos.band_u_sw[lv,ba] = max(0.0,atmos.radout.flux_up[idx])
@@ -661,12 +661,13 @@ module energy
     - `conduct::Bool`                   include conductive heat transport
     - `convect_sf::Float64`             scale factor applied to convection fluxes
     - `latent_sf::Float64`              scale factor applied to phase change fluxes
-    - `calc_cf::Bool=false`             calculate LW contribution function?
+    - `calc_cf::Bool`                   calculate LW contribution function?
+    - `reset_vmrs::Bool`                reset VMRs to dry values before radtrans and MLT
     """
     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)
+                          calc_cf::Bool=false, reset_vmrs::Bool=true)
 
         # Reset fluxes
         reset_fluxes!(atmos)
@@ -689,9 +690,11 @@ module energy
             @turbo @. atmos.flux_tot += atmos.flux_l
 
             # Restore mixing ratios
-            for g in atmos.gas_names
-                @turbo @. atmos.gas_vmr[g] = atmos.gas_ovmr[g]
-            end
+            if reset_vmrs
+                for g in atmos.gas_names
+                    @turbo @. atmos.gas_vmr[g] = atmos.gas_ovmr[g]
+                end
+            end 
         end
         # Calculate layer properties
         atmosphere.calc_layer_props!(atmos)

src/plotting.jl

@@ -100,6 +100,44 @@ module plotting
         return plt
     end
 
+    """
+    Plot the height vs pressure profile.
+    """
+    function plot_height(atmos::atmosphere.Atmos_t, fname::String;
+                                dpi::Int=250,
+                                size_x::Int=500, size_y::Int=400,
+                                incl_magma::Bool=false,
+                                title::String="")
+
+        ylims  = (1e-5*atmos.pl[1]/1.5, 1e-5*atmos.pl[end]*1.5)
+        yticks = 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))
+
+        # Create plot
+        plt = plot(ylims=ylims, yticks=yticks, legend=:outertopright,
+                        dpi=dpi, size=(size_x,size_y))
+
+        # Plot surface
+        scatter!(plt, [0.0], [atmos.pl[end]*1e-5], color="brown3", label=L"P_s")
+
+        # Plot cell-centres and cell-edges
+        scatter!(plt, atmos.z*1e-3,  atmos.p*1e-5,  msa=0.0, msw=0, ms=1.2, shape=:diamond, label="Centres")
+        scatter!(plt, atmos.zl*1e-3, atmos.pl*1e-5, msa=0.0, msw=0, ms=1.2, shape=:diamond, label="Edges")
+
+        # Decorate
+        xlabel!(plt, "Height [km]")
+        ylabel!(plt, "Pressure [bar]")
+        yflip!(plt)
+        yaxis!(plt, yscale=:log10)
+        if !isempty(title)
+            title!(plt, title)
+        end
+
+        if !isempty(fname)
+            savefig(plt, fname)
+        end
+        return plt
+    end
+
     """
     Plot the cloud mass mixing ratio and area fraction.
     """

src/solver.jl

@@ -76,8 +76,9 @@ module solver
     - `sol_type::Int`                   solution type, 1: tmp_surf | 2: skin | 3: flux_int | 4: tgt_olr
     - `chem_type::Int`                  chemistry type (see wiki)
     - `convect::Bool`                   include convection
-    - `sens_heat::Bool`                 include sensible heating
+    - `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)
     - `dx_max::Float64`                 maximum step size [K]
     - `max_steps::Int`                  maximum number of solver steps
     - `max_runtime::Float64`            maximum runtime in wall-clock seconds
@@ -103,8 +104,8 @@ module solver
     function solve_energy!(atmos::atmosphere.Atmos_t;
                             sol_type::Int=1,
                             chem_type::Int=0,
-                            convect::Bool=true, sens_heat::Bool=false,
-                            conduct::Bool=false, latent::Bool=false,
+                            convect::Bool=true, sens_heat::Bool=true,
+                            conduct::Bool=true, latent::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,

src/spectrum.jl

@@ -108,8 +108,18 @@ module spectrum
             @error "Minimum wavelength is too small"
             return false
         end
+        if wl[2] < wl[1]
+            @error "Stellar wavelength array must be strictly ascending"
+            return false
+        end
         clamp!(fl, 1.0e-45, 1.0e+45)  # Clamp values
 
+        # Ensure that wl array has no duplicates
+        # https://discourse.julialang.org/t/unique-indices-method-similar-to-matlab/34446/3
+        unique_mask::Array{Int} = unique(z -> wl[z], 1:length(wl))
+        wl = wl[unique_mask]
+        fl = fl[unique_mask]
+
         # Bin data to required number of bins...
 
         # Log data first