Skip to content

Radiation

frierson_atmospheric_heating(ds, albedo=0)

Returns the atmospheric radiative heating rate from the surface and top of atmosphere energy fluxes. A negative value indicates that the atmosphere is cooling.

This takes into account any radiation that is absorbed by the atmosphere on its way down from space to the surface, as specified through tau_equator and the amount reflected at the surface through the albedo.

In Isca, there is no absorption of the shortwave radiation as it moves back up through the atmosphere to space after being reflected at the surface.

Parameters:

Name Type Description Default
ds Dataset

Dataset for particular experiment, must contain:

  • swdn_toa - Incident shortwave radiation at the top of atmosphere. This is saved by Isca if the variable swdn_toa in the two_stream module is specified in the diagnostic table.
  • swdn_sfc - Net shortwave radiation absorbed at the surface i.e. incident - reflected. This is the negative of net upward shortwave radiation at the surface. This is saved by Isca if the variable swdn_sfc in the two_stream module is specified in the diagnostic table.
  • lwup_sfc - Upward longwave flux at the surface. This is saved by Isca if the variable lwdn_sfc in the two_stream module is specified in the diagnostic table.
  • lwdn_sfc - Downward longwave flux at the surface. This is saved by Isca if the variable lwdn_sfc in the two_stream module is specified in the diagnostic table.
  • olr - Outgoing longwave radiation at the top of atmosphere. This is saved by Isca if the variable lwdn_sfc in the two_stream module is specified in the diagnostic table.
 required
albedo float

Fraction of incident shortwave radiation reflected by the surface. It is specified through the option albedo_value in the mixed_layer_nml namelist.

 0

Returns: Atmospheric radiative heating rate in \(W/m^2\).

Source code in isca_tools/utils/radiation.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def frierson_atmospheric_heating(ds: Dataset, albedo: float = 0) -> xr.DataArray:
    """
    Returns the atmospheric radiative heating rate from the surface and top of atmosphere energy fluxes. A negative
    value indicates that the atmosphere is cooling.

    This takes into account any radiation that is absorbed by the atmosphere on its way down from space to the surface,
    as specified through `tau_equator` and the amount reflected at the surface through the `albedo`.

    In *Isca*, there is no absorption of the shortwave radiation as it moves back up through the atmosphere to space
    after being reflected at the surface.

    Args:
        ds: Dataset for particular experiment, must contain:

            * `swdn_toa` - Incident shortwave radiation at the top of atmosphere.
                This is saved by *Isca* if the variable `swdn_toa` in the `two_stream` module is specified in the
                diagnostic table.
            * `swdn_sfc` - Net shortwave radiation absorbed at the surface i.e. incident - reflected.
                This is the negative of net upward shortwave radiation at the surface.
                This is saved by *Isca* if the variable `swdn_sfc` in the `two_stream` module is specified in the
                diagnostic table.
            * `lwup_sfc` - Upward longwave flux at the surface.
                This is saved by *Isca* if the variable `lwdn_sfc` in the `two_stream` module is specified in the
                diagnostic table.
            * `lwdn_sfc` - Downward longwave flux at the surface.
                This is saved by *Isca* if the variable `lwdn_sfc` in the `two_stream` module is specified in the
                diagnostic table.
            * `olr` - Outgoing longwave radiation at the top of atmosphere.
                This is saved by *Isca* if the variable `lwdn_sfc` in the `two_stream` module is specified in the
                diagnostic table.
        albedo: Fraction of incident shortwave radiation reflected by the surface.
            It is specified through the option `albedo_value` in the `mixed_layer_nml` namelist.
    Returns:
        Atmospheric radiative heating rate in $W/m^2$.

    """
    # #This is the full method of doing it, but we can do it simpler without the sw absorption stuff
    # swup_net_sfc = -ds.swdn_sfc
    # # need to account for SW absorbed by atmosphere and reflected at surface
    # swup_net_toa = -frierson_net_toa_sw_dwn(ds.swdn_toa, ds.ps, albedo, tau_equator, tau_lat_var, pressure_exponent,
    #                                         ref_pressure)
    # flux_surf = swup_net_sfc + ds.lwup_sfc - ds.lwdn_sfc
    # flux_toa = swup_net_toa + ds.olr
    # return ds.swdn_toa - ds.swdn_sfc / (1 - albedo) + ds.lwup_sfc - ds.lwdn_sfc - ds.olr

    swup_toa = albedo / (1 - albedo) * ds.swdn_sfc
    flux_toa = swup_toa - ds.swdn_toa + ds.olr
    flux_surf = -ds.swdn_sfc + ds.lwup_sfc - ds.lwdn_sfc
    return flux_surf - flux_toa

frierson_net_toa_sw_dwn(insolation, surface_pressure, albedo=0, tau_equator=0, tau_lat_var=0, pressure_exponent=4, ref_pressure=101325)

Function to calculate the net downward shortwave radiation at the top of atmosphere for the Frierson and Byrne radiation schemes.

This takes into account any radiation that is absorbed by the atmosphere on its way down from space to the surface, as specified through tau_equator and the amount reflected at the surface through the albedo.

In Isca, there is no absorption of the shortwave radiation as it moves back up through the atmosphere to space after being reflected at the surface.

All default values are the default values in Isca if the option is not specified in the relavent namelist.

Parameters:

Name Type Description Default
insolation DataArray

Incident shortwave radiation at the top of atmosphere with dimensions of latitude, longitude and time. This is saved by Isca if the variable swdn_toa in the two_stream module is specified in the diagnostic table.

 required
surface_pressure DataArray

Surface pressure in Pa with dimensions of latitude, longitude and time, \(p_s\). This is saved by Isca if the variable ps in the dynamics module is specified in the diagnostic table.

 required
albedo float

Fraction of incident shortwave radiation reflected by the surface. It is specified through the option albedo_value in the mixed_layer_nml namelist.

 0
tau_equator float

Surface optical depth at the equator, \(\tau_{0e}^*\). It is specified through the option atm_abs in the two_stream_gray_rad_nml namelist.

 0
tau_lat_var float

Variation of optical depth with latitude, \(\Delta \tau^*\). It is specified through the option sw_diff in the two_stream_gray_rad_nml namelist.

 0
pressure_exponent float

Determines the variation of optical depth with pressure, \(\kappa^*\). It is specified through the option solar_exponent in the two_stream_gray_rad_nml namelist.

 4
ref_pressure float

Reference pressure used by Isca in Pa. It is specified through the option pstd_mks in the constants_nml namelist.

 101325

Returns:

Type Description
DataArray

Net downward shortwave radiation at the top of atmosphere with dimensions of latitude, longitude and time.
Source code in isca_tools/utils/radiation.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def frierson_net_toa_sw_dwn(insolation: xr.DataArray, surface_pressure: xr.DataArray, albedo: float = 0,
                            tau_equator: float = 0, tau_lat_var: float = 0, pressure_exponent: float = 4,
                            ref_pressure: float = 101325) -> xr.DataArray:
    """
    Function to calculate the net downward shortwave radiation at the top of atmosphere for the *Frierson*
    and *Byrne* radiation schemes.

    This takes into account any radiation that is absorbed by the atmosphere on its way down from space to the surface,
    as specified through `tau_equator` and the amount reflected at the surface through the `albedo`.

    In *Isca*, there is no absorption of the shortwave radiation as it moves back up through the atmosphere to space
    after being reflected at the surface.

    All default values are the default values in *Isca* if the option is not specified in the relavent `namelist`.

    Args:
        insolation: Incident shortwave radiation at the top of atmosphere
            with dimensions of latitude, longitude and time.
            This is saved by *Isca* if the variable `swdn_toa` in the `two_stream` module is specified in the
            diagnostic table.
        surface_pressure: Surface pressure in *Pa* with dimensions of latitude, longitude and time, $p_s$.
            This is saved by *Isca* if the variable `ps` in the `dynamics` module is specified in the diagnostic table.
        albedo: Fraction of incident shortwave radiation reflected by the surface.
            It is specified through the option `albedo_value` in the `mixed_layer_nml` namelist.
        tau_equator: Surface optical depth at the equator, $\\tau_{0e}^*$.
            It is specified through the option `atm_abs` in the `two_stream_gray_rad_nml` namelist.
        tau_lat_var: Variation of optical depth with latitude, $\Delta \\tau^*$.
            It is specified through the option `sw_diff` in the `two_stream_gray_rad_nml` namelist.
        pressure_exponent: Determines the variation of optical depth with pressure, $\\kappa^*$.
            It is specified through the option `solar_exponent` in the `two_stream_gray_rad_nml` namelist.
        ref_pressure: Reference pressure used by Isca in *Pa*.
            It is specified through the option `pstd_mks` in the `constants_nml` namelist.

    Returns:
        Net downward shortwave radiation at the top of atmosphere with dimensions of latitude, longitude and time.
    """
    tau = frierson_sw_optical_depth(surface_pressure, tau_equator, tau_lat_var, pressure_exponent, ref_pressure)
    return insolation * (1 - albedo * np.exp(-tau))

frierson_sw_optical_depth(surface_pressure, tau_equator=0, tau_lat_var=0, pressure_exponent=4, ref_pressure=101325)

Function to calculate shortwave surface optical depth, \(\tau_s\), as a function of latitude, \(\phi\), as performed in Isca for the Frierson and Byrne radiation schemes:

\[ \tau_s(\phi) = (1-\Delta \tau^* \sin^2 \phi)\tau_{0e}^*\big(\frac{p_s}{p_0}\big)^{\kappa^*} \]

All default values are the default values in Isca if the option is not specified in the relavent namelist.

Parameters:

Name Type Description Default
surface_pressure DataArray

Surface pressure in Pa with dimensions of latitude, longitude and time, \(p_s\). This is saved by Isca if the variable ps in the dynamics module is specified in the diagnostic table.

 required
tau_equator float

Surface optical depth at the equator, \(\tau_{0e}^*\). It is specified through the option atm_abs in the two_stream_gray_rad_nml namelist.

 0
tau_lat_var float

Variation of optical depth with latitude, \(\Delta \tau^*\). It is specified through the option sw_diff in the two_stream_gray_rad_nml namelist.

 0
pressure_exponent float

Determines the variation of optical depth with pressure, \(\kappa^*\). It is specified through the option solar_exponent in the two_stream_gray_rad_nml namelist.

 4
ref_pressure float

Reference pressure used by Isca in Pa. It is specified through the option pstd_mks in the constants_nml namelist.

 101325

Returns:

Type Description
DataArray

Shortwave surface optical depth with dimensions of latitude, longitude and time.
Source code in isca_tools/utils/radiation.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def frierson_sw_optical_depth(surface_pressure: xr.DataArray, tau_equator: float = 0, tau_lat_var: float = 0,
                              pressure_exponent: float = 4, ref_pressure: float = 101325) -> xr.DataArray:
    """
    Function to calculate shortwave surface optical depth, $\\tau_s$, as a function of latitude, $\\phi$, as performed in
    [Isca](https://execlim.github.io/Isca/modules/two_stream_gray_rad.html#frierson-byrne-schemes) for the *Frierson*
    and *Byrne* radiation schemes:

    $$
    \\tau_s(\\phi) = (1-\Delta \\tau^* \sin^2 \\phi)\\tau_{0e}^*\\big(\\frac{p_s}{p_0}\\big)^{\\kappa^*}
    $$

    All default values are the default values in *Isca* if the option is not specified in the relavent `namelist`.

    Args:
        surface_pressure: Surface pressure in *Pa* with dimensions of latitude, longitude and time, $p_s$.
            This is saved by *Isca* if the variable `ps` in the `dynamics` module is specified in the diagnostic table.
        tau_equator: Surface optical depth at the equator, $\\tau_{0e}^*$.
            It is specified through the option `atm_abs` in the `two_stream_gray_rad_nml` namelist.
        tau_lat_var: Variation of optical depth with latitude, $\Delta \\tau^*$.
            It is specified through the option `sw_diff` in the `two_stream_gray_rad_nml` namelist.
        pressure_exponent: Determines the variation of optical depth with pressure, $\\kappa^*$.
            It is specified through the option `solar_exponent` in the `two_stream_gray_rad_nml` namelist.
        ref_pressure: Reference pressure used by Isca in *Pa*.
            It is specified through the option `pstd_mks` in the `constants_nml` namelist.

    Returns:
        Shortwave surface optical depth with dimensions of latitude, longitude and time.
    """
    tau_surface = (1 - tau_lat_var * np.sin(np.deg2rad(surface_pressure.lat)) ** 2) * tau_equator
    return tau_surface * (surface_pressure / ref_pressure) ** pressure_exponent

get_frierson_sw_abs(atm_abs=None, p_surf=None, sw_diff=0, solar_exponent=4, p_ref=101325, swdn_sfc=None, swdn_toa=None, albedo=0)

Calculates the fraction of incoming shortwave radiation absorbed by the atmosphere.

Calculates atmospheric shortwave absorption either from supplied downward shortwave fluxes at the surface and top of atmosphere, or from the Frierson shortwave optical-depth parameterization. The flux-based method is used only when both swdn_sfc and swdn_toa are provided.

Parameters:

Name Type Description Default
atm_abs Optional[float]

Atmospheric shortwave absorption parameter used by the optical-depth parameterization. Ignored when both shortwave flux inputs are provided. This is the name in Isca namelist.

 None
p_surf Optional[DataArray]

Surface pressure. Used by the optical-depth parameterization; ignored when both shortwave flux inputs are provided.

 None
sw_diff float

Additive shortwave optical-depth parameter. This is the name in Isca namelist.

 0
solar_exponent float

Exponent controlling the pressure dependence of shortwave optical depth. This is the name in Isca namelist.

 4
p_ref float

Reference pressure for the optical-depth parameterization, in Pa.

 101325
swdn_sfc Optional[DataArray]

Downward shortwave flux at the surface.

 None
swdn_toa Optional[DataArray]

Downward shortwave flux at the top of the atmosphere.

 None
albedo float

Surface albedo used to infer atmospheric absorption from the fluxes.

 0

Returns:

Type Description
DataArray

Fraction of incoming shortwave radiation absorbed by the atmosphere.
DataArray

The returned xarray.DataArray has the dimensions and coordinates of
DataArray

the supplied fluxes or surface pressure.
Notes

With supplied fluxes, atmospheric absorption is

\[ f_{\mathrm{abs}} = 1 - \frac{F_{\mathrm{sfc}}} {F_{\mathrm{toa}}(1 - \alpha)}. \]

Otherwise, it is calculated from the shortwave optical depth as

\[ f_{\mathrm{abs}} = 1 - \exp(-\tau_{\mathrm{sw}}), \]

where $\(\tau_{\mathrm{sw}}\)$ is calculated by frierson_sw_optical_depth.

Source code in isca_tools/utils/radiation.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def get_frierson_sw_abs(atm_abs: Optional[float] = None, p_surf: Optional[xr.DataArray] = None,
                        sw_diff: float = 0, solar_exponent: float = 4, p_ref: float = 101325,
                        swdn_sfc: Optional[xr.DataArray] = None, swdn_toa: Optional[xr.DataArray] = None,
                        albedo: float = 0) -> xr.DataArray:
    r"""Calculates the fraction of incoming shortwave radiation absorbed by the atmosphere.

    Calculates atmospheric shortwave absorption either from supplied downward
    shortwave fluxes at the surface and top of atmosphere, or from the
    Frierson shortwave optical-depth parameterization. The flux-based method
    is used only when both ``swdn_sfc`` and ``swdn_toa`` are provided.

    Args:
        atm_abs: Atmospheric shortwave absorption parameter used by the
            optical-depth parameterization. Ignored when both shortwave flux
            inputs are provided. This is the name in Isca namelist.
        p_surf: Surface pressure. Used by the optical-depth parameterization;
            ignored when both shortwave flux inputs are provided.
        sw_diff: Additive shortwave optical-depth parameter. This is the name in Isca namelist.
        solar_exponent: Exponent controlling the pressure dependence of
            shortwave optical depth. This is the name in Isca namelist.
        p_ref: Reference pressure for the optical-depth parameterization, in
            Pa.
        swdn_sfc: Downward shortwave flux at the surface.
        swdn_toa: Downward shortwave flux at the top of the atmosphere.
        albedo: Surface albedo used to infer atmospheric absorption from the
            fluxes.

    Returns:
        Fraction of incoming shortwave radiation absorbed by the atmosphere.
        The returned `xarray.DataArray` has the dimensions and coordinates of
        the supplied fluxes or surface pressure.

    Notes:
        With supplied fluxes, atmospheric absorption is

        $$
        f_{\mathrm{abs}} =
        1 - \frac{F_{\mathrm{sfc}}}
        {F_{\mathrm{toa}}(1 - \alpha)}.
        $$

        Otherwise, it is calculated from the shortwave optical depth as

        $$
        f_{\mathrm{abs}} = 1 - \exp(-\tau_{\mathrm{sw}}),
        $$

        where $$\tau_{\mathrm{sw}}$$ is calculated by
        `frierson_sw_optical_depth`.
    """
    if (swdn_sfc is not None) and (swdn_toa is not None):
        # I think this method may have problems with diurnal cycle if using daily average data
        sw_abs = 1 - swdn_sfc / swdn_toa / (1 - albedo)
    else:
        odp_sw = frierson_sw_optical_depth(p_surf, atm_abs, sw_diff, solar_exponent, p_ref)
        sw_abs = 1 - np.exp(-odp_sw)  # fraction of sw absorbed
    return sw_abs

get_heat_capacity(c_p, density, layer_depth, return_depth=False)

Given heat capacity un units of \(JK^{-1}kg^{-1}\), this returns heat capacity in units of \(JK^{-1}m^{-2}\).

Parameters:

Name Type Description Default
c_p float

Specific heat at constant pressure.
Units: \(JK^{-1}kg^{-1}\)

 required
density float

Density of substance (usually air or water).
Units: \(kgm^{-3}\)

 required
layer_depth float

Depth of layer or heat capacity if return_depth=True
Units: \(m\)

 required
return_depth bool

If True, will return mixed layer depth in \(m\) given heat capacity in \(JK^{-1}m^{-2}\).

 False

Returns:

Type Description
float

Heat capacity in units of \(JK^{-1}m^{-2}\) or mixed layer depth in \(m\) if return_depth=True.
Source code in isca_tools/utils/radiation.py 
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
def get_heat_capacity(c_p: float, density: float, layer_depth: float, return_depth: bool = False) -> float:
    """
    Given heat capacity un units of $JK^{-1}kg^{-1}$, this returns heat capacity in units of $JK^{-1}m^{-2}$.

    Args:
        c_p: Specific heat at constant pressure.</br>
            Units: $JK^{-1}kg^{-1}$
        density: Density of substance (usually air or water).</br>
            Units: $kgm^{-3}$
        layer_depth: Depth of layer or heat capacity if `return_depth=True`</br>
            Units: $m$
        return_depth: If `True`, will return mixed layer depth in $m$ given heat capacity in $JK^{-1}m^{-2}$.

    Returns:
        Heat capacity in units of $JK^{-1}m^{-2}$ or mixed layer depth in $m$ if `return_depth=True`.
    """
    if return_depth:
        return layer_depth / (c_p * density)
    else:
        return c_p * density * layer_depth

get_sw_abs_amp(swdn_sfc, swdn_toa, time, albedo)

Compute the atmospheric shortwave absorption from annual-harmonic amplitudes.

The atmospheric absorption fraction is inferred from the ratio between the amplitude of the first annual harmonic of downward shortwave radiation at the surface and at the top of the atmosphere:

\(f_\mathrm{abs} = 1 - \frac{A_\mathrm{sfc}} {A_\mathrm{TOA}(1 - \alpha)}\)

where \(A_\mathrm{sfc}\) and \(A_\mathrm{TOA}\) are the positive amplitudes of the annual harmonic in surface and top-of-atmosphere downward shortwave radiation, respectively, and \(\alpha\) is the surface albedo.

Parameters:

Name Type Description Default
swdn_sfc ndarray

Downward shortwave radiation at the surface, sampled at time.

 required
swdn_toa ndarray

Downward shortwave radiation at the top of the atmosphere, sampled at time.

 required
time ndarray

Time coordinate corresponding to the radiation data. Its units and sampling must be compatible with :func:get_fourier_coef.

 required
albedo float

Surface albedo, expressed as a fraction between 0 and 1.

 required

Returns:

Type Description
float

Fraction of the incoming annual-harmonic shortwave-radiation amplitude
float

absorbed by the atmosphere.
Notes

This diagnostic concerns the seasonal-cycle amplitude rather than the time-mean shortwave absorption. It assumes that amp_toa and 1 - albedo are non-zero.

Source code in isca_tools/utils/radiation.py 
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def get_sw_abs_amp(swdn_sfc: np.ndarray, swdn_toa: np.ndarray, time: np.ndarray, albedo: float) -> float:
    r"""Compute the atmospheric shortwave absorption from annual-harmonic amplitudes.

    The atmospheric absorption fraction is inferred from the ratio between the
    amplitude of the first annual harmonic of downward shortwave radiation at
    the surface and at the top of the atmosphere:

    $f_\mathrm{abs} = 1 - \frac{A_\mathrm{sfc}}
    {A_\mathrm{TOA}(1 - \alpha)}$

    where $A_\mathrm{sfc}$ and $A_\mathrm{TOA}$ are the positive amplitudes of
    the annual harmonic in surface and top-of-atmosphere downward shortwave
    radiation, respectively, and $\alpha$ is the surface albedo.

    Args:
        swdn_sfc:
            Downward shortwave radiation at the surface, sampled at `time`.
        swdn_toa:
            Downward shortwave radiation at the top of the atmosphere, sampled at
            `time`.
        time:
            Time coordinate corresponding to the radiation data. Its units and
            sampling must be compatible with :func:`get_fourier_coef`.
        albedo:
            Surface albedo, expressed as a fraction between 0 and 1.

    Returns:
            Fraction of the incoming annual-harmonic shortwave-radiation amplitude
            absorbed by the atmosphere.

    Notes:
        This diagnostic concerns the seasonal-cycle amplitude rather than the
        time-mean shortwave absorption. It assumes that `amp_toa` and
        `1 - albedo` are non-zero.
    """
    amp_sfc = get_fourier_coef(time, swdn_sfc, 1, pos_amp=True)[0]
    amp_toa = get_fourier_coef(time, swdn_toa, 1, pos_amp=True)[0]
    sw_abs = 1 - amp_sfc / amp_toa / (1 - albedo)
    return sw_abs

opd_lw_gray(lat, pressure=None, kappa=1, tau_eq=6, tau_pole=1.5, pressure_ref=10 ** 5, frac_linear=0.1, k_exponent=4)

Returns the longwave optical depth used in the Frierson Isca rad_scheme.

If pressure not provided, will return surface value. Otherwise, will return value at give pressure.

Parameters:

Name Type Description Default
lat ndarray

float [n_lat]
Latitude in degrees.

 required
pressure Optional[float]

Pressure in Pa.

 None
kappa float

Frierson optical depth scaling parameter.
opd in two_stream_gray_rad_nml namelist.

 1
tau_eq float

Surface longwave optical depth at equator.
ir_tau_eq in two_stream_gray_rad_nml namelist.

 6
tau_pole float

Surface longwave optical depth at pole.
ir_tau_pole in two_stream_gray_rad_nml namelist.

 1.5
pressure_ref float

Reference pressure in Pa.

 10 ** 5
frac_linear float

Determines partitioning between linear term and \(p^k\) term.
linear_tau in two_stream_gray_rad_nml namelist.

 0.1
k_exponent float

Pressure exponent.
wv_exponent in two_stream_gray_rad_nml namelist.

 4

Returns:

Name Type Description
opd ndarray

float [n_lat]
Longwave optical depth
Source code in isca_tools/utils/radiation.py 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
def opd_lw_gray(lat: np.ndarray, pressure: Optional[float] = None,
                kappa: float = 1, tau_eq: float = 6, tau_pole: float = 1.5,
                pressure_ref: float = 10 ** 5, frac_linear: float = 0.1, k_exponent: float = 4) -> np.ndarray:
    """
    Returns the longwave optical depth used in the
    [Frierson](https://execlim.github.io/Isca/modules/two_stream_gray_rad.html#frierson-byrne-schemes)
    Isca `rad_scheme`.

    If `pressure` not provided, will return surface value. Otherwise, will return value at give `pressure`.

    Args:
        lat: `float [n_lat]`</br>
            Latitude in degrees.
        pressure: Pressure in Pa.
        kappa: Frierson optical depth scaling parameter.</br>
            `opd` in `two_stream_gray_rad_nml` namelist.
        tau_eq: Surface longwave optical depth at equator.</br>
            `ir_tau_eq` in `two_stream_gray_rad_nml` namelist.
        tau_pole: Surface longwave optical depth at pole.</br>
            `ir_tau_pole` in `two_stream_gray_rad_nml` namelist.
        pressure_ref: Reference pressure in Pa.
        frac_linear: Determines partitioning between linear term and $p^k$ term.</br>
            `linear_tau` in `two_stream_gray_rad_nml` namelist.
        k_exponent: Pressure exponent.</br>
            `wv_exponent` in `two_stream_gray_rad_nml` namelist.

    Returns:
        opd: `float [n_lat]`</br>
            Longwave optical depth
    """
    opd_surf = kappa * (tau_eq + (tau_pole - tau_eq) * np.sin(np.deg2rad(lat)) ** 2)
    if pressure is None:
        return opd_surf
    else:
        pressure_factor = frac_linear * (pressure / pressure_ref) + (1 - frac_linear) * (
                pressure / pressure_ref) ** k_exponent
        return opd_surf * pressure_factor