You can download the tar.gz archive of the source code here.
This archive corresponds to the version used by E. Marcq et al. in their 2017 paper available here.
This code was used to produce the results discussed in Marcq et al. (2017) paper. It simulates a 1D atmospheric column of H2O and CO2 above a solid surface, and computes a low-resolution spectrum of the thermal emission at the top of the atmosphere. For more information about the physics of the model, please refer to the above cited paper:
- E. Marcq, A. Salvador, H. Massol and A. Davaille, Thermal radiation of magma ocean planets using a 1D radiative-convective model of H2O-CO2 atmospheres, JGR (2017)
It makes extensive use of the provided and publicly available DISORT radiative transfer code:
- K. Stamnes et al., Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl. Opt. 27 (1988)
A working Fortran 90/95 compiler is required. Our code has only been tested using the following Fortran compilers and OS:
- pgfortran 13.10-0 64-bit under Scientific Linux release 6.8 (x86_64)
- gfortran 5.4.0 under Ubuntu 16.04.2 LTS (x86_64)
Although not strictly necessary, it is strongly recommended to have a working f2py installation, available in the NumPy library for Python 2.7+. This enables to use the model through any Python v2.7+ environment with Numpy installed.
Source code tree
sub_radconv.f95 | ---make_profiles.f95 | | | ---steam.f90 | ---make_tauk.f95 | | | ---k_correle.f95 | | | ---make_tau_grey.f95 | ---make_tau_grey.f95 | ---disort_wrapper.f95 | ---(all provided DISORT subroutines)
- g-grid for k-correlated data
- T-grid for k-correlated data
- p-grid for k-correlated data
- Q-grid for k-correlated data
- H2O-H2O continuum LUT from MTCKD
- CO2-CO2 continuum LUT from B. Bézard
- k-correlated coefficients computed by KSPECTRUM
- directory required for diagnostic Fortran outputs (created empty by unpacking the *.tar.gz file, so it should be OK).
- untar the archive file
- compile the source files
- using the
compil_python.shscript (if you have a working
f2pyinstallation) to produce a Python library named
radconv1d.sothat you can import from Python — you may want to change the path to your python installation in this script.
- using the
compil_fortran.shscript for a Fortran-only installation. This will produce a
radconv1dexecutable file. The main program is just a wrapper which provides the main subroutine the required inputs (see below for a description of these inputs).
- using the
Using the model with Python (recommended)
Example 1 (main subroutine)
from radconv1d import radconv Ttop_out, flux_up, P0_v, P0_c, P0_o, spectrum = radconv(Tsurf, Ttop_in, Psurf_v, Psurf_c, Psurf_o, Ttop_flag, cloud_flag, cont_flag, grey_flag, GP_flag, spher_flag, alb_clear, alb_cloud, const_sol, g0, Rt, Nz)
- Surface temperature in K
- Starting mesospheric temperature in K. 200 is a recommanded value.
- Total surface pressure of H2O in Pa (including liquid H2O if any)
- Total surface pressure of CO2 in Pa
- Total surface pressure of N2 in Pa
- if == 1, Ttop will not be altered. Else, the model will try to adjust Ttop so that the TOA divergence of OLR is null. This results in unrealistically cold mesospheres, and is kept only in anticipation of stellar heating modeling!
- if == 0, clouds won't be radiatively active. If == 1, Earth-like clouds are radiatively assumed. If == 2, Venus-like clouds are radiatively assumed (not realistic for H2O-rich atmospheres)
- if == 0, no continuum opacity is included. If == 1, only CO2-CO2 opacity is included. If == 2, only H2O-H2O opacity is included. If == 3, H2O-H2O and CO2-CO2 opacities are included (recommended).
- if == 0, line opacities are computed according to k-correlated data (slow, but more realistic). Else, line opacities will be averaged over the whole spectrum (much faster, but less accurate for the total integrated OLR, and plainly wrong for the OLR spectrum).
- if ==1, H2O steam is treated as an ideal gas. Else, its EOS is given by the included steam tables.
- future flag for sphericity corrections. Should be set equal to 0.
- Prescribed albedo without clouds. Preliminary SW calculations yield about 0.2 around a G-star.
- Prescribed albedo with clouds. Preliminary SW calculations yield about 0.7 around a G-star.
- Averaged stellar constant in W/m². == 1/4 of its value at the substellar point.
- Surface gravity in m/s²
- Planetary radius in m
- Number of computational layers. Typically around 200.
- Final mesospheric temperature
in K (differs from
Ttop_flag != 1).
- Total OLR in W/m²
- Gaseous surface pressure of H2O in Pa (excluding surface liquid H2O if any).
- Gaseous surface pressure of CO2 in Pa (== Psurf_c in this version)
- Gaseous surface pressure of N2 in Pa (== Psurf_o in this version)
- Spectrum of OLR in W/m²/cm-1 in each of the 36 bands. Boundary wavenumbers are given below:
sigmas = np.array([0., 40., 160., 280., 380., 500., 582., 600., 720., 752., 800., 900., 1000., 1200., 1350., 1450., 1550., 1650., 1750., 1850., 1950., 2050., 2200., 2500., 2800., 3200., 3600., 4000., 4400., 4800., 5900.,6000., 6500., 8000., 8300., 9300., 10100.]) # in inv-cm
Example 2 (
from radconv1d import make_profiles Z, T, Pv, Pc, Po, rho_v, rho_c, rho_o, ibottom, itop = make_profiles(GP_flag, g0, Rt, Tsurf, Psurf_v, Psurf_c, Psurf_o, Ttop, Nz)
Input parametersSee above.
(Nz+1)-long 1D array of boundary altitudes for computational layers in m.
(Nz+1)-long 1D array of boundary temperatures for computational layers in K.
(Nz+1)-long 1D array of boundary H2O pressures for computational layers in Pa.
(Nz+1)-long 1D array of boundary CO2 pressures for computational layers in Pa.
(Nz+1)-long 1D array of boundary N2 pressures for computational layers in Pa.
(Nz)-long array of H2O densities within computational layers in kg/m3.
(Nz)-long array of CO2 densities within computational layers in kg/m3.
(Nz)-long array of N2 densities within computational layers in kg/m3.
- index of the lowermost computational layer in the moist adiabat layer
- index of the uppermost computational layer in the moist adiabat layer
Using the model with Fortran only
Starting from the provided
Fortran wrapper, the user is free to provide required inputs to
the main subroutine named
RADCONV, and process its
ouputs afterwards. If no spectral output is required but merely
vertical profiles in composition, density and temperature, you
can call only
For the physical meaning of the inputs and outputs of the aforementioned subroutines, please refer to the Python section hereabove.
If you find an error in the code or wish to submit a feature request, please contact the author by e-mail.
Dr Emmanuel Marcq, Associate Professor, LATMOS -
Univ. Paris-Saclay Tel: +33 1 80 28 52 83
This wode has been released under Creative Common License BY-SA-NC
BY: Licensees may copy, distribute, display and perform the work
and make derivative works and remixes based on it only if they
give the author or licensor the credits (attribution) in the
"Author: E. Marcq/LATMOS-IPSL-UPSAY (email@example.com)"
SA: Licensees may distribute derivative works only under a license identical ("not more restrictive") to the license that governs the original work.
NC: Licensees may copy, distribute, display, and perform the work and make derivative works and remixes based on it only for non-commercial purposes.
- To K. Stamnes et al. for writing and making publicly available the DISORT code
- To J. Leconte and R. Wordsworth for their help during the testing and development of this code.