Merge pull request #77 from nichollsh/hn/contfunc

Contribution function updates

CITATION.cff

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

Project.toml

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

src/AGNI.jl

@@ -443,7 +443,7 @@ module AGNI
         @info prt_req[1:end-2]
 
         # Write initial state
-        dump.write_ptz(atmos, joinpath(atmos.OUT_DIR,"ptz_ini.csv"))
+        dump.write_ptz(atmos, joinpath(atmos.OUT_DIR,"ptz_initial.csv"))
 
         # Do chemistry on initial composition
         if chem_type in [1,2,3]
@@ -525,8 +525,8 @@ module AGNI
 
         # Write arrays
         @info "Writing results"
-        dump.write_ptz(atmos,      joinpath(atmos.OUT_DIR,"ptz.csv"))
-        dump.write_fluxes(atmos,  joinpath(atmos.OUT_DIR,"fl.csv"))
+        # dump.write_ptz(atmos,     joinpath(atmos.OUT_DIR,"ptz.csv"))
+        # dump.write_fluxes(atmos,  joinpath(atmos.OUT_DIR,"fl.csv"))
         dump.write_ncdf(atmos,    joinpath(atmos.OUT_DIR,"atm.nc"))
 
         # Save plots
@@ -537,7 +537,6 @@ module AGNI
         plt_ani && plotting.animate(atmos)
         plt_vmr && plotting.plot_vmr(atmos,       joinpath(atmos.OUT_DIR,"plot_vmrs.png"), size_x=600)
         plt_cff && plotting.plot_contfunc1(atmos, joinpath(atmos.OUT_DIR,"plot_contfunc1.png"))
-        plt_cff && plotting.plot_contfunc2(atmos, joinpath(atmos.OUT_DIR,"plot_contfunc2.png"))
         plt_tmp && plotting.plot_pt(atmos,        joinpath(atmos.OUT_DIR,"plot_ptprofile.png"), incl_magma=(sol_type==2))
         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"))

src/atmosphere.jl

@@ -13,7 +13,7 @@ module atmosphere
     using LinearAlgebra
     using Logging
     using LoopVectorization
-    # import Statistics
+    import Statistics
     import PCHIPInterpolation:Interpolator
     import DelimitedFiles:readdlm
 
@@ -209,7 +209,7 @@ module atmosphere
         fastchem_work::String           # Path to fastchem working directory
 
         # Observing properties
-        transspec_p::Float64            # (INPUT PARAMETER) pressure level probed in transmission [Pa]
+        transspec_p::Float64            # Pressure level probed in transmission [Pa]
         transspec_r::Float64            # planet radius probed in transmission [m]
         transspec_m::Float64            # mass [kg] enclosed by transspec_r
         transspec_rho::Float64          # bulk density [kg m-3] implied by r and m
@@ -333,7 +333,7 @@ module atmosphere
         @info  "Setting-up a new atmosphere struct"
 
         # Code versions
-        atmos.AGNI_VERSION = "0.11.1"
+        atmos.AGNI_VERSION = "0.11.2"
         atmos.SOCRATES_VERSION = readchomp(joinpath(ENV["RAD_DIR"],"version"))
         @debug "AGNI VERSION = "*atmos.AGNI_VERSION
         @debug "Using SOCRATES at $(ENV["RAD_DIR"])"
@@ -717,10 +717,64 @@ module atmosphere
         return true
     end # end function setup
 
+
+    """
+    **Estimate photosphere.**
+
+    Estimates the location of the photosphere by finding the median of the contribution
+    function in each band (0.2 um to 150 um), and then finding the pressure level at which
+    these median values are maximised.
+
+    Arguments:
+        - `atmos::Atmos_t`          the atmosphere struct instance to be used.
+
+    Returns:
+        - `p_ref::Float64`          pressure level of photosphere [Pa]
+    """
+    function estimate_photosphere!(atmos::atmosphere.Atmos_t)::Float64
+
+        # Params
+        wl_min::Float64  = 0.2 * 1e-6 # 200 nm
+        wl_max::Float64  = 150 * 1e-6 # 150 um
+        p_min::Float64   = 10.0 # 1e-4 bar
+
+        # tracking
+        cff_max::Float64 = 0.0
+        cff_try::Float64 = 0.0
+        atmos.transspec_p = p_min
+
+        # get band indices
+        wl_imin = findmin(abs.(atmos.bands_cen .- wl_min))[2]
+        wl_imax = findmin(abs.(atmos.bands_cen .- wl_max))[2]
+
+        # reversed?
+        if wl_imin > wl_imax
+            wl_imin, wl_imax = wl_imax, wl_imin
+        end
+
+        # loop over levels
+        for i in 1:atmos.nlev_c
+            if atmos.p[i] < p_min
+                continue
+            end
+
+            # maximum contfunc in this band
+            cff_try = Statistics.median(atmos.contfunc_band[i,wl_imin:wl_imax])
+
+            # is this more than the existing maximum?
+            if cff_try > cff_max
+                cff_max = cff_try
+                atmos.transspec_p = atmos.p[i]
+            end
+        end
+
+        return atmos.transspec_p
+    end
+
     """
     **Calculate observed radius and bulk density.**
 
-    This is done at the layer probed in transmission, which is set at a fixed pressure.
+    This is done at the layer probed in transmission, which is set to a fixed pressure.
 
     Arguments:
         - `atmos::Atmos_t`          the atmosphere struct instance to be used.
@@ -729,15 +783,14 @@ module atmosphere
     Returns:
         - `transspec_rho::Float64`  the bulk density observed in transmission
     """
-    function calc_observed_rho!(atmos::atmosphere.Atmos_t, p_ref::Float64=100.0)::Float64
+    function calc_observed_rho!(atmos::atmosphere.Atmos_t)::Float64
 
-        # transspec_p::Float64            # (INPUT) level probed in transmission [Pa]
         # transspec_r::Float64            # planet radius probed in transmission [m]
         # transspec_m::Float64            # mass [kg] of atmosphere + interior
         # transspec_rho::Float64          # bulk density [kg m-3] implied by r and m
 
         # Store reference pressure in atmos struct
-        atmos.transspec_p = p_ref
+        estimate_photosphere!(atmos)
 
         # get the observed height
         idx::Int = findmin(abs.(atmos.p .- atmos.transspec_p))[2]

src/dump.jl

@@ -125,6 +125,7 @@ module dump
             # ----------------------
             # Scalar quantities
             #    Create variables
+            var_max_cff_p = defVar(ds, "max_cff_p",     Float64, (), attrib = OrderedDict("units" => "Pa"))     # Pressure level at which CFF is max [Pa]
             var_tmp_surf =  defVar(ds, "tmp_surf",      Float64, (), attrib = OrderedDict("units" => "K"))      # Surface brightness temperature [K]
             var_flux_int =  defVar(ds, "flux_int",      Float64, (), attrib = OrderedDict("units" => "W m-2"))  # Internal flux [W m-2]
             var_inst =      defVar(ds, "instellation",  Float64, (), attrib = OrderedDict("units" => "W m-2"))  # Solar flux at TOA
@@ -149,6 +150,7 @@ module dump
             var_starfile =  defVar(ds, "starfile"      ,String, ())     # Path to star file when read
 
             #     Store data
+            var_max_cff_p[1] =  atmos.transspec_p
             var_tmp_surf[1] =   atmos.tmp_surf
             var_flux_int[1] =   atmos.flux_int
             var_inst[1] =       atmos.instellation