Better interpolation.

Project.toml

@@ -1,7 +1,7 @@
 name = "AGNI"
 uuid = "ede838c1-9ec3-4ebe-8ae8-da4091b3f21c"
 authors = ["Harrison Nicholls <[email protected]>"]
-version = "1.1.0"
+version = "1.1.1"
 
 [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": "1.1.0"
+  "version": "1.1.1"
 }

res/config/structure.toml

@@ -0,0 +1,70 @@
+# AGNI configuration file
+title = "Hot and dry"
+
+[planet]
+    tmp_surf        = 2000.0
+    instellation    = 44000.0
+    albedo_b        = 0.18
+    s0_fact         = 0.6652
+    zenith_angle    = 60.0
+    surface_material= "res/surface_albedos/lunar_marebasalt.dat"
+    albedo_s        = 0.0
+    radius          = 6.37e6
+    gravity         = 9.81
+    flux_int        = 0.0
+    turb_coeff      = 1.0e-2
+    wind_speed      = 10.0
+
+[files]
+    input_sf        = "res/spectral_files/Frostflow/256/Frostflow.sf"
+    input_star      = "res/stellar_spectra/sun.txt"
+    output_dir      = "out/"
+
+[composition]
+    p_surf          = 2700.0
+    p_top           = 1e-5
+    vmr_dict        = { H2O = 1.0}
+    vmr_path        = ""
+    include_all     = false
+    chemistry       = 0
+    condensates     = []
+
+
+[execution]
+    clean_output    = true
+    verbosity       = 1
+    max_steps       = 200
+    max_runtime     = 400
+    num_levels      = 50
+    continua        = true
+    rayleigh        = true
+    cloud           = false
+    aerosol         = false
+    overlap_method  = "ee"
+    real_gas        = false
+    thermo_funct    = true
+    gravity_funct   = true
+    sensible_heat   = true
+    latent_heat     = false
+    convection      = true
+    rainout         = false
+    solution_type   = 1
+    solver          = ""
+    dx_max          = 400.0
+    initial_state   = ["dry","sat","H2O"]
+    linesearch      = 0
+    easy_start      = false
+    converge_atol   = 1.0e-2
+    converge_rtol   = 1.0e-4
+
+
+[plots]
+    at_runtime      = true
+    temperature     = true
+    fluxes          = true
+    contribution    = true
+    emission        = true
+    albedo          = true
+    mixing_ratios   = true
+    animate         = true
+    height          = true

src/atmosphere.jl

@@ -14,7 +14,7 @@ module atmosphere
     using Logging
     using LoopVectorization
     import Statistics
-    import Interpolations: interpolate, Gridded, Linear, Flat, extrapolate, Extrapolation
+    import Interpolations: interpolate, Gridded, Linear, Flat, Line, extrapolate, Extrapolation
     import DelimitedFiles:readdlm
 
     # Local files
@@ -336,7 +336,7 @@ module atmosphere
         @info  "Setting-up a new atmosphere struct"
 
         # Code versions
-        atmos.AGNI_VERSION = "1.1.0"
+        atmos.AGNI_VERSION = "1.1.1"
         atmos.SOCRATES_VERSION = readchomp(joinpath(ENV["RAD_DIR"],"version"))
         @debug "AGNI VERSION = "*atmos.AGNI_VERSION
         @debug "Using SOCRATES at $(ENV["RAD_DIR"])"
@@ -1036,14 +1036,14 @@ module atmosphere
         # Logarithmically-spaced levels above
         atmos.pl[1:end-1] .= collect(Float64, range( start=atmos.pl[1],
                                                       stop=atmos.pl[end-1],
-                                                      length=atmos.nlev_l-1))
+                                                            length=atmos.nlev_l-1))
 
         # Shrink near-surface layers by stretching all layers above
         p_fact::Float64 = 0.6
         p_mid::Float64 = atmos.pl[end-1]*p_fact + atmos.pl[end-2]*(1.0-p_fact)
         atmos.pl[1:end-2] .= collect(Float64, range( start=atmos.pl[1],
-                                                     stop=p_mid,
-                                                     length=atmos.nlev_l-2))
+                                                            stop=p_mid,
+                                                            length=atmos.nlev_l-2))
 
         # Set cell-centres at midpoint of cell-edges
         atmos.p[1:end] .= 0.5 .* (atmos.pl[1:end-1] .+ atmos.pl[2:end])
@@ -2122,17 +2122,10 @@ module atmosphere
     function set_tmpl_from_tmp!(atmos::atmosphere.Atmos_t)
 
         # Interpolate temperature to bulk cell-edge values (log-linear)
-        itp = extrapolate(interpolate((log10.(atmos.p),),atmos.tmp,Gridded(Linear())),Flat())
-        atmos.tmpl[2:end-1] .= itp.(log10.(atmos.pl[2:end-1]))
-
-        # Extrapolate top edge temperature (log-linear)
-        grad_dt::Float64 = atmos.tmp[1] - atmos.tmp[2]
-        grad_dp::Float64 = log10(atmos.p[1]/atmos.p[2])
-        atmos.tmpl[1] = atmos.tmp[1] + grad_dt/grad_dp * log10(atmos.pl[1]/atmos.p[1])
-
-        # Set bottom edge to bottom cell-centre value
-        # This is fine because the bottom cell is very small (in pressure space)
-        atmos.tmpl[end]=atmos.tmp[end]
+        itp = extrapolate( interpolate( (log10.(atmos.p[:]),),
+                                        atmos.tmp,Gridded(Linear())
+                                        ),Line())
+        atmos.tmpl[:] .= itp.(log10.(atmos.pl))
 
         # Clamp
         clamp!(atmos.tmpl, atmos.tmp_floor, atmos.tmp_ceiling)

src/solver.jl

@@ -407,10 +407,10 @@ module solver
 
         # Allocate initial guess for the x array, as well as a,b arrays
         # Array storage structure:
-        #   in the sol_type=2 case
+        #   in the sol_type>=2 cases
         #       1:end-1 => cell centre temperatures
         #       end     => bottom cell edge temperature
-        #   other cases
+        #   sol_type == 1 case
         #       1:end => cell centre temperatures
         x_ini    = zeros(Float64, arr_len)
         for i in 1:atmos.nlev_c