In [1]:
import os

os.environ['ARTS_DATA_PATH'] = '/home/pc2943/arts-cat-data/'
In [2]:
import matplotlib.pyplot as plt
from typhon import plots

import konrad


plots.styles.use('typhon')

/home/pc2943/.local/lib/python3.9/site-packages/typhon/nonlte/rtc/__init__.py:9: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def FOSC(tau, Sb, Sm, Ib):
In [3]:
import numpy as np
import pyarts
from typhon import physics as phys
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.pyplot as plt
from pathlib import Path
import pandas as pd

import copy
import seaborn as sns
import colorcet as cc

from scipy.constants import speed_of_light as c

import xarray as xr
import pyarts.workspace

import re

Documentation¶

Quadrature¶

The quadrature scheme included as a radiation option is as described in Czarnecki, P., Polvani, L., & Pincus, R. (2023). Sparse, empirically optimized quadrature for broadband spectral integration. Journal of Advances in Modeling Earth Systems, 15(10), e2023MS003819.

There are two schemes, one for the shortwave and one for the longwave, included. These are both trained to parameterize variations in present-day amounts of water vapor, ozone, and temperature throughout the atmosphere as top-of-the-atmosphere carbon dioxide forcing in clear skies. This scheme is not yet suitable for perturbations in other gases and has not been tested with clouds.

The default scheme remains RRTMG. The quadrature scheme can be used by calling the following inside konrad.RCE():

radiation = konrad.radiation.ARTS(                                          
    arts_kwargs = {                                                           
        'quadrature':True,                                                    
        'quadrature_filename_lw': '/path/to/lw_quadrature_configuration.h5',  
        'lookup_filename_lw': '/path/to/lw_lookup_table.xml',                 
        'quadrature_filename_sw': '/path/to/sw_quadrature_configuration.h5',  
        'lookup_filename_sw': '/path/to/sw_lookup_table.xml',                 
    }                                                                         
)

If 'quadrature':True, the quadrature scheme is used. Otherwise, line-by-line ARTS calculations are performed across about 60,000 frequencies in the longwave and the shortwave. quardature_filename_lw and quadrature_filename_sw are netCDF files that contain at least the following variables: an array $S$ with the wavenumbers (in units of cm$^{-1}$) of the quadrature points in ascending order, and an array $W$ of the corresponding weights, in the corresponding order. If undefined or None, the default configuration, stored in konrad/konrad/radiation/data, will be used. lookup_filename_lw and lookup_filename_sw are by default the gas absorption lookup tables corresponding to the default quadrature configuration. If set as None, line-by-line calculations will be performed by ARTS for the given quadrature points.

ARTS¶

Line-by-line ARTS calculations can also be performed in the shortwave and the longwave. This is achieved by calling the following inside konrad.RCE():

radiation = konrad.radiation.ARTS()

There are a few requirements to running ARTS, both line-by-line and with the quadrature scheme:

  • pyARTS version 2.5.11 was used to implement ARTS/quadrature in konrad.
  • Lookup tables to the ARTS line-by-line calculation must be provided. This can be accomplished by running the following in the command line before activating the python environment:
    export KONRAD_LOOKUP_TABLE_LW='/path/to/abs_lookup_lw.xml'
    export KONRAD_LOOKUP_TABLE_SW='/path/to/abs_lookup_sw.xml'
    An example longwave lookup table can be found at https://doi.org/10.5281/zenodo.3885410.
  • The path to the ARTS catalogue must be provided. This can be accomplished by running the following line in Python before importing ARTS or konrad:
    import os
    os.environ['ARTS_DATA_PATH'] = '/path/to/arts-cat-data/'
    The konrad lookup tables can be similarly included.
    The catalogue must include:
    • line data (in a folder called lines)
    • cross-section induced absorption (CIA) data (with the path cia-xml/hitran2011/hitran_cia2012_adapted.xml.gz)
    • cross-section data (in a folder called xsec)
    • MT_CKD 4.0 water vapor contiuum data (with the path model/mt_ckd_4.0/H2O.xml)
    • the 2008 July solar spectrum for the shortwave (with the path arts-xml-data/star/Sun/solar_spectrum_July_2008_konrad_scale.xml).
      Most of these can be found in the arts-cat-data repository (https://arts.mi.uni-hamburg.de/svn/rt/arts-cat-data/branches/). In particular, ozone cross-sections must be provided for the shortwave to work correctly.

Making Lookup Tables¶

In order to make lookup tables for either a quadrature scheme configuration or line-by-line ARTS calculations, one can use built-in functions for either the longwave or the shortwave by running the following code: 

# create an instance of the ARTS class
arts_class = konrad.radiation.arts._ARTS()

# create sw lookup table
konrad.radiation.arts._ARTS.calc_lookup_table_sw(
    arts_class, 
    filename = 'path/to/lookup_filename.xml', 
    t_min = min_temperature,                     # default 150 K
    t_max = max_temperature,                     # default 350 K
    p_step = pressure_discretization_step_size,  # default 0.5
    wavenumber = wavenumber                      # default 2^15 wavenumbers linearly spaced btw. 2,000 - 50,000 cm^{-1}
)

# or for lw
konrad.radiation.arts._ARTS.calc_lookup_table_lw(
    arts_class, 
    filename = 'path/to/lookup_filename.xml', 
    f_num = number_of_wavenumbers,               # default 2^15 linearly spaced btw. 10 - 3,250 cm^{-1}
    wavenumber = wavenumber                      # caution! Input in units m^{-1}
)

Examples¶

Below are several examples demonstrating the working of the quadrature scheme as compared to the default RRTMG radiation code and the ARTS line-by-line code. Most come from the konrad documentation notebook (https://atmtools.github.io/konrad/getting_started.html).

Present-Day Radiative Convective Equilibrium¶

Here, we run several RCE simulations with fixed surface temperatures that span the range of present-day conditions. The temperature profile, net flux profile, and heating rates that result from the final converged state of the model are plotted for the three different radiation codes - the default RRTMG, the line-by-line ARTS code, and the Quadrature scheme. The Quadrature scheme tracks ARTS well for these present-day conditions.

In [4]:
temps = np.array([220, 250, 280, 310])
for T in temps:
    
    # define the quadrature atmosphere
    plev, phlev = konrad.utils.get_pressure_grids(1000e2, 1, 128)
    quad_atmosphere = konrad.atmosphere.Atmosphere(phlev)
    
    # run the RCE model with quadrature as the radiation scheme
    quad = konrad.RCE(
        quad_atmosphere,
        surface=konrad.surface.FixedTemperature(temperature=T),  # Run with a fixed surface temperature.
        timestep='12h',  # Set timestep in model time.
        max_duration='100d',  # Set maximum runtime.
        radiation = konrad.radiation.ARTS(arts_kwargs = {'quadrature':True})
    )
    
    quad.run()
    
    # To run the same ARTS simulation, one can use this code. 
    # It will take some time.
    '''
    # define the ARTS atmosphere
    plev, phlev = konrad.utils.get_pressure_grids(1000e2, 1, 128)
    arts_atmosphere = konrad.atmosphere.Atmosphere(phlev)
    
    # run konrad with ARTS as the radiation calculation
    arts = konrad.RCE(
        arts_atmosphere,
        surface=konrad.surface.FixedTemperature(temperature=T),  # Run with a fixed surface temperature.
        timestep='12h',  # Set timestep in model time.
        max_duration='100d',  # Set maximum runtime.
        radiation = konrad.radiation.ARTS()
    )
    
    
    arts.run()
    
    # save results to file    
    nc = konrad.netcdf.NetcdfHandler('arts_output' + str(T) + '_sw.nc', arts)  # create output file
    nc.write()  # write (append) current RCE state to file'''

    arts_atmosphere = konrad.atmosphere.Atmosphere.from_netcdf('arts_output_' + str(T) +'_quadrature.nc')
    arts_surface = konrad.surface.FixedTemperature.from_netcdf('arts_output_' + str(T) +'_quadrature.nc')
    arts_class = konrad.radiation.arts.ARTS()
    konrad.radiation.arts.ARTS.calc_radiation(arts_class, arts_atmosphere, arts_surface, konrad.cloud.ClearSky)
    konrad.radiation.arts.ARTS.update_heatingrates(arts_class, arts_atmosphere, arts_surface)
    arts_radiation = arts_class.data_vars
       
    
    # define the rrtmg atmosphere    
    plev, phlev = konrad.utils.get_pressure_grids(1000e2, 1, 128)
    rrtmg_atmosphere = konrad.atmosphere.Atmosphere(phlev)

    # run the default RCE model
    rrtmg = konrad.RCE(
        rrtmg_atmosphere,
        surface=konrad.surface.FixedTemperature(temperature=T),  # Run with a fixed surface temperature.
        timestep='12h',  # Set timestep in model time.
        max_duration='100d',  # Set maximum runtime.
    )
    
    rrtmg.run()
        
    # define net fluxes
    lw_net_rrtmg = rrtmg.radiation.data_vars['lw_flxd'][1][0] - rrtmg.radiation.data_vars['lw_flxu'][1][0]
    lw_net_arts = arts_radiation['lw_flxd'][1][0] - arts_radiation['lw_flxu'][1][0]
    lw_net_quad = quad.radiation.data_vars['lw_flxd'][1][0] - quad.radiation.data_vars['lw_flxu'][1][0]

    sw_net_rrtmg = rrtmg.radiation.data_vars["sw_flxd"][1][0] - rrtmg.radiation.data_vars["sw_flxu"][1][0]
    sw_net_arts = arts_radiation["sw_flxd"][1][0] - arts_radiation["sw_flxu"][1][0]
    sw_net_quad = quad.radiation.data_vars["sw_flxd"][1][0] - quad.radiation.data_vars["sw_flxu"][1][0]

    ### PLOTS
    fig = plt.figure()
    fig.suptitle('T$_s$ = ' + str(T) + ' K', fontsize = 16)

    ### Temperature    
    ax = plt.subplot2grid((2, 2), (0, 0), rowspan = 2)
    ax.semilogy(rrtmg_atmosphere['T'][-1, :], rrtmg_atmosphere['plev']/100, label = 'rrtmg', color = 'teal')
    ax.semilogy(arts_atmosphere['T'][-1, :], arts_atmosphere['plev']/100, label = 'arts', color = 'maroon')
    ax.semilogy(quad_atmosphere['T'][-1, :], quad_atmosphere['plev']/100, label = 'quadrature', color = 'gold')

    ax.set_ylim([1100, 0.01])
    #ax.legend(frameon = False)
    ax.set_xlabel(r"$T$ (K)")
    ax.set_ylabel("$p$ (hPa)")

    ### Fluxes
    
    ax = plt.subplot2grid((2, 2), (0, 1))
    ax.semilogy(lw_net_rrtmg, rrtmg_atmosphere['phlev']/100, label = 'rrtmg', color = 'teal')
    ax.semilogy(lw_net_arts, arts_atmosphere['phlev']/100, label = 'arts', color = 'maroon')
    ax.semilogy(lw_net_quad, quad_atmosphere['phlev']/100, label = 'quadrature', color = 'gold')
    
    ax.semilogy(sw_net_rrtmg, rrtmg_atmosphere['phlev']/100, '--', label = '_nolegend_', color = 'teal')
    ax.semilogy(sw_net_arts, arts_atmosphere['phlev']/100, '--', label = '_nolegend_', color = 'maroon')
    ax.semilogy(sw_net_quad, quad_atmosphere['phlev']/100, '--', label = '_nolegend_', color = 'gold')

    ax.set_ylim([1100, 0.01])

    ax.plot(0, 1200, '--', label = 'shortwave',  color = 'black')
    ax.plot(0, 1200, label = 'longwave',  color = 'black')

    ax.legend(frameon = False)
    
    ax.set_ylabel('pressure (hPa)')
    ax.set_xlabel('net radiative flux ($\sf W/m^2$)')

    ### Heating Rates
    
    ax = plt.subplot2grid((2, 2), (1, 1))
    ax.axvline(0, color='black', linewidth=0.8)
    ax.semilogy(rrtmg.radiation.data_vars["lw_htngrt"][1][0], rrtmg_atmosphere['plev']/100, label = 'rrtmg', color = 'teal')
    ax.semilogy(arts_radiation["lw_htngrt"][1][0], arts_atmosphere['plev']/100, label = 'arts', color = 'maroon')
    ax.semilogy(quad.radiation.data_vars["lw_htngrt"][1][0], quad_atmosphere['plev']/100, label = 'quadrature', color = 'gold')

    ax.semilogy(rrtmg.radiation.data_vars["sw_htngrt"][1][0], rrtmg_atmosphere['plev']/100, '--', label = '_nolegend_', color = 'teal')
    ax.semilogy(arts_radiation["sw_htngrt"][1][0], arts_atmosphere['plev']/100, '--', label = '_nolegend_', color = 'maroon')
    ax.semilogy(quad.radiation.data_vars["sw_htngrt"][1][0], quad_atmosphere['plev']/100, '--', label = '_nolegend_', color = 'gold')

    ax.set_ylim([1100, 0.001])
    

    plt.xlabel("heating rate (K/day)")
    plt.ylabel("pressure (hPa)")
    plt.show()

Instantaneous TOA Forcing by CO$_2$¶

We calculate instantaneous top-of-the-atmosphere forcing by CO$_2$ as described by the konrad guide (https://atmtools.github.io/konrad/forcing.html). The atmosphere is allowed to equilibriate for present-day CO$_2$ concentrations. Then the CO$_2$ concentration is uniformly perturbed throughout the atmosphere and fluxes and heating rates are re-calculated. Forcing profiles and temperatures through the atmosphere are reported for models run with the three radiation schemes. The quadrature is trained on 8xCO$_2$ and can interpolate to different forcing states, so we show 2x, 4x, and 8x CO$_2$.

First, we allow the present-day atmosphere to equilibriate.

In [4]:
phlev = konrad.utils.get_quadratic_pgrid(1000e2, 10, 128)
atmosphere = konrad.atmosphere.Atmosphere(phlev)
atmosphere["CO2"][:] = 348e-6 # Set reference CO2 concentration

# Calculate reference OLR.
spinup = konrad.RCE(atmosphere, 
                    timestep='24h', 
                    max_duration='150d',     
                    radiation = konrad.radiation.ARTS(arts_kwargs = {'quadrature':True})
)
spinup.run()
In [5]:
phlev = konrad.utils.get_quadratic_pgrid(1000e2, 10, 128)
atmosphere_rrtmg = konrad.atmosphere.Atmosphere(phlev)
atmosphere_rrtmg["CO2"][:] = 348e-6 # Set reference CO2 concentration

# Calculate reference OLR.
spinup_rrtmg = konrad.RCE(atmosphere_rrtmg, timestep='24h', max_duration='150d')
spinup_rrtmg.run()

2xCO$_2$¶

Next, we double the CO$_2$ concentration, update the heating rates, and calculate the instantaneous forcing.

In [6]:
olr_ref = spinup.radiation["lw_flxu"][-1]
# Calculate OLR at perturbed atmospheric state.
atmosphere["CO2"][:] *= 2  # double the CO2 concentration
spinup.radiation.update_heatingrates(atmosphere, spinup.surface)

instant_forcing = -(spinup.radiation["lw_flxu"][-1] - olr_ref)
In [7]:
olr_ref_rrtmg = spinup_rrtmg.radiation["lw_flxu"][-1]

# Calculate OLR at perturbed atmospheric state.
atmosphere_rrtmg["CO2"][:] *= 2  # double the CO2 concentration
spinup_rrtmg.radiation.update_heatingrates(atmosphere_rrtmg, spinup_rrtmg.surface)

instant_forcing_rrtmg = -(spinup_rrtmg.radiation["lw_flxu"][-1] - olr_ref_rrtmg)

The ARTS calculation is loaded from file, since the convergence takes some time.

In [8]:
arts_atmosphere = konrad.atmosphere.Atmosphere.from_netcdf('arts_output_co2.nc')
arts_surface = konrad.surface.FixedTemperature.from_netcdf('arts_output_co2.nc')
arts_class = konrad.radiation.arts.ARTS()
konrad.radiation.arts.ARTS.calc_radiation(arts_class, arts_atmosphere, arts_surface, konrad.cloud.ClearSky)
konrad.radiation.arts.ARTS.update_heatingrates(arts_class, arts_atmosphere, arts_surface)
arts_radiation = arts_class.data_vars
In [9]:
olr_ref_arts = arts_radiation["lw_flxu"][-1]

# Calculate OLR at perturbed atmospheric state.
arts_atmosphere["CO2"][:] *= 2  # double the CO2 concentration
konrad.radiation.arts.ARTS.update_heatingrates(arts_class, arts_atmosphere, arts_surface)

instant_forcing_arts = -(arts_radiation["lw_flxu"][-1] - olr_ref_arts)

Finally, we plot the results.

In [10]:
fig, (ax0, ax1) = plt.subplots(ncols=2, sharey=True)
fig.suptitle('2xCO$_2$', fontsize = 16)

plots.profile_p_log(spinup_rrtmg.atmosphere["phlev"], instant_forcing_rrtmg, ax=ax0, color = 'teal')
plots.profile_p_log(arts_atmosphere["phlev"], instant_forcing_arts[0], ax = ax0, color = 'maroon')
plots.profile_p_log(spinup.atmosphere["phlev"], instant_forcing, ax=ax0, color = 'gold')

ax0.set_xticks([0, np.round(instant_forcing_arts[0][-1], 2), np.round(instant_forcing_arts[0].max(), 2)])
ax0.grid(axis="x")
ax0.set_xlim(0, 4.5)
ax0.set_xlabel(r"$\Delta F$ (W m$^{-2}$)")
ax0.set_ylabel("$p$ (hPa)")
ax0.set_ylim(bottom=1000e2)
ax0.legend(['rrtmg', 'arts', 'quadrature'])

plots.profile_p_log(spinup_rrtmg.atmosphere["plev"], spinup_rrtmg.atmosphere["T"][-1], ax=ax1, color = 'teal')
plots.profile_p_log(arts_atmosphere["plev"], arts_atmosphere["T"][-1], ax = ax1, color = 'maroon')
plots.profile_p_log(spinup.atmosphere["plev"], spinup.atmosphere["T"][-1], ax=ax1, color = 'gold')

ax1.set_xticks([196, 267, 288])
ax1.grid(axis="x")
ax1.set_xlim(180, 290)
ax1.set_xlabel(r"$T$ (K)")

print(f"Instantanoues TOA forcing calculated by Quadrature: {instant_forcing[-1]:.2f} W/m^2")
print(f"Instantanoues TOA forcing calculated by RRTMG: {instant_forcing_rrtmg[-1]:.2f} W/m^2")
print(f"Instantanoues TOA forcing calculated by ARTS: {instant_forcing_arts[0][-1]:.2f} W/m^2")

Instantanoues TOA forcing calculated by Quadrature: 2.58 W/m^2
Instantanoues TOA forcing calculated by RRTMG: 2.57 W/m^2
Instantanoues TOA forcing calculated by ARTS: 2.52 W/m^2

4xCO$_2$¶

The previous CO$_2$ concentration is again doubled for all models.

In [11]:
# Calculate OLR at perturbed atmospheric state.
atmosphere_rrtmg["CO2"][:] *= 2  # double the CO2 concentration
spinup_rrtmg.radiation.update_heatingrates(atmosphere_rrtmg, spinup_rrtmg.surface)

instant_forcing_rrtmg = -(spinup_rrtmg.radiation["lw_flxu"][-1] - olr_ref_rrtmg)

atmosphere["CO2"][:] *= 2  # double the CO2 concentration
spinup.radiation.update_heatingrates(atmosphere, spinup.surface)

instant_forcing = -(spinup.radiation["lw_flxu"][-1] - olr_ref)

# Calculate OLR at perturbed atmospheric state.
arts_atmosphere["CO2"][:] *= 2  # double the CO2 concentration
konrad.radiation.arts.ARTS.update_heatingrates(arts_class, arts_atmosphere, arts_surface)

instant_forcing_arts = -(arts_radiation["lw_flxu"][-1] - olr_ref_arts)

fig, (ax0, ax1) = plt.subplots(ncols=2, sharey=True)
fig.suptitle('4xCO$_2$', fontsize = 16)

plots.profile_p_log(spinup_rrtmg.atmosphere["phlev"], instant_forcing_rrtmg, ax=ax0, color = 'teal')
plots.profile_p_log(arts_atmosphere["phlev"], instant_forcing_arts[0], ax = ax0, color = 'maroon')
plots.profile_p_log(spinup.atmosphere["phlev"], instant_forcing, ax=ax0, color = 'gold')

ax0.set_xticks([0, np.round(instant_forcing_arts[0][-1], 2), np.round(instant_forcing_arts[0].max(), 2)])
ax0.grid(axis="x")
#ax0.set_xlim(0, 4.5)
ax0.set_xlabel(r"$\Delta F$ (W m$^{-2}$)")
ax0.set_ylabel("$p$ (hPa)")
ax0.set_ylim(bottom=1000e2)
ax0.legend(['rrtmg', 'arts', 'quadrature'])

plots.profile_p_log(spinup_rrtmg.atmosphere["plev"], spinup_rrtmg.atmosphere["T"][-1], ax=ax1, color = 'teal')
plots.profile_p_log(arts_atmosphere["plev"], arts_atmosphere["T"][-1], ax = ax1, color = 'maroon')
plots.profile_p_log(spinup.atmosphere["plev"], spinup.atmosphere["T"][-1], ax=ax1, color = 'gold')

ax1.set_xticks([196, 267, 288])
ax1.grid(axis="x")
ax1.set_xlim(180, 290)
ax1.set_xlabel(r"$T$ (K)")

print(f"Instantanoues TOA forcing calculated by Quadrature: {instant_forcing[-1]:.2f} W/m^2")
print(f"Instantanoues TOA forcing calculated by RRTMG: {instant_forcing_rrtmg[-1]:.2f} W/m^2")
print(f"Instantanoues TOA forcing calculated by ARTS: {instant_forcing_arts[0][-1]:.2f} W/m^2")

Instantanoues TOA forcing calculated by Quadrature: 5.06 W/m^2
Instantanoues TOA forcing calculated by RRTMG: 5.21 W/m^2
Instantanoues TOA forcing calculated by ARTS: 5.05 W/m^2

8xCO$_2$¶

We double CO$_2$ again and update the radiation calculation.

In [12]:
# Calculate OLR at perturbed atmospheric state.
atmosphere_rrtmg["CO2"][:] *= 2  # double the CO2 concentration
spinup_rrtmg.radiation.update_heatingrates(atmosphere_rrtmg, spinup_rrtmg.surface)

instant_forcing_rrtmg = -(spinup_rrtmg.radiation["lw_flxu"][-1] - olr_ref_rrtmg)

atmosphere["CO2"][:] *= 2  # double the CO2 concentration
spinup.radiation.update_heatingrates(atmosphere, spinup.surface)

instant_forcing = -(spinup.radiation["lw_flxu"][-1] - olr_ref)

# Calculate OLR at perturbed atmospheric state.
arts_atmosphere["CO2"][:] *= 2  # double the CO2 concentration
konrad.radiation.arts.ARTS.update_heatingrates(arts_class, arts_atmosphere, arts_surface)

instant_forcing_arts = -(arts_radiation["lw_flxu"][-1] - olr_ref_arts)


fig, (ax0, ax1) = plt.subplots(ncols=2, sharey=True)
fig.suptitle('8xCO$_2$', fontsize = 16)

plots.profile_p_log(spinup_rrtmg.atmosphere["phlev"], instant_forcing_rrtmg, ax=ax0, color = 'teal')
plots.profile_p_log(arts_atmosphere["phlev"], instant_forcing_arts[0], ax = ax0, color = 'maroon')
plots.profile_p_log(spinup.atmosphere["phlev"], instant_forcing, ax=ax0, color = 'gold')

ax0.set_xticks([0, np.round(instant_forcing_arts[0][-1], 2), np.round(instant_forcing_arts[0].max(), 2)])
ax0.grid(axis="x")
#ax0.set_xlim(0, 4.5)
ax0.set_xlabel(r"$\Delta F$ (W m$^{-2}$)")
ax0.set_ylabel("$p$ (hPa)")
ax0.set_ylim(bottom=1000e2)
ax0.legend(['rrtmg', 'arts', 'quadrature'])

plots.profile_p_log(spinup_rrtmg.atmosphere["plev"], spinup_rrtmg.atmosphere["T"][-1], ax=ax1, color = 'teal')
plots.profile_p_log(arts_atmosphere["plev"], arts_atmosphere["T"][-1], ax = ax1, color = 'maroon')
plots.profile_p_log(spinup.atmosphere["plev"], spinup.atmosphere["T"][-1], ax=ax1, color = 'gold')

ax1.set_xticks([196, 267, 288])
ax1.grid(axis="x")
ax1.set_xlim(180, 290)
ax1.set_xlabel(r"$T$ (K)")

print(f"Instantanoues TOA forcing calculated by Quadrature: {instant_forcing[-1]:.2f} W/m^2")
print(f"Instantanoues TOA forcing calculated by RRTMG: {instant_forcing_rrtmg[-1]:.2f} W/m^2")
print(f"Instantanoues TOA forcing calculated by ARTS: {instant_forcing_arts[0][-1]:.2f} W/m^2")

Instantanoues TOA forcing calculated by Quadrature: 7.46 W/m^2
Instantanoues TOA forcing calculated by RRTMG: 7.82 W/m^2
Instantanoues TOA forcing calculated by ARTS: 7.60 W/m^2
In [ ]:
In [ ]: