r"""
Classes to provide insolation as input for other CLIMLAB processes.
Options include
- ``climlab.radiation.P2Insolation`` (idealized 2nd Legendre polynomial form)
- ``climlab.radiation.FixedInsolation`` (generic steady-state insolation)
- ``climlab.radiation.AnnualMeanInsolation`` (steady-state annual-mean insolation computed from orbital parameters and latitude)
- ``climlab.radiation.DailyInsolation`` (time-varying daily-mean insolation computed from orbital parameters, latitude and time of year)
- ``climlab.radiation.InstantInsolation`` (time-varying instantaneous insolation computed from orbital parameters, latitude, longitude, and time of year)
All are subclasses of ``climlab.process.DiagnosticProcess``
and do not add any tendencies to any state variables.
At least three diagnostics are provided:
- ``insolation``, the incoming solar radiation in :math:`\frac{\textrm{W}}{\textrm{m}^2}`
- ``coszen``, cosine of the solar zenith angle (dimensionless)
- ``irradiance_factor``, ratio of current total irradiance to its annual average, i.e. solar constant (dimensionless)
"""
import numpy as np
from climlab.process.diagnostic import DiagnosticProcess
from climlab.domain.field import Field, to_latlon
from climlab.utils.legendre import P2
from climlab import constants as const
from climlab.solar.insolation import daily_insolation_factors, instant_insolation_factors, \
annual_insolation, dates_to_day_index
# REVISE TO MAKE ALL OF THESE CALLABLE WITH NO ARGUMENTS.
# SET SOME SENSIBLE DEFAULTS FOR DOMAINS
# the diagnostic self.insolation is set with correct dimensions at
# creation time. After that, make sure to always modify it though
# self.insolation[:] = ...
# so that links to the insolation in other processes will work
# REALLY NEED TO CHANGE THE WAY DIAGNOSTIC PROCESSES ARE CREATED
# CAN'T RELY ON NAMES OF DOMAINS
# should be easy to pass a state variable object or something with the right shape
# and have the process gracefully set the correct dimensions
class _Insolation(DiagnosticProcess):
r"""A private parent class for insolation processes.
Calling compute() will update self.insolation with current values.
**Initialization parameters** \n
An instance of ``_Insolation`` is initialized with the following
arguments *(for detailed information see Object attributes below)*:
:param float S0: solar constant
- unit: :math:`\frac{\textrm{W}}{\textrm{m}^2}`
- default value: ``1365.2``
**Object attributes**
Additional to the parent class :class:`~climlab.process.diagnostic.DiagnosticProcess`
following object attributes are generated and updated during initialization:
:ivar array insolation: the array is initialized with zeros of the size of
``self.domains['sfc']`` or ``self.domains['default']``.
:ivar array coszen: cosine of the solar zenith angle
:ivar float S0: initialized with given argument ``S0``
:ivar dict diagnostics: key ``'insolation'`` initialized with value:
:class:`~climlab.domain.field.Field` of zeros
in size of ``self.domains['sfc']`` or
``self.domains['default']``
:ivar Field insolation: the subprocess attribute ``self.insolation`` is
created with correct dimensions
.. note::
``self.insolation`` should always be modified with
``self.insolation[:] = ...`` so that links to the insolation in other
processes will work.
"""
# parameter S0 is now stored using a python property
# can be changed through self.S0 = newvalue
# which will also update the parameter dictionary
# CAUTION: changing self.param['S0'] will not work!
def __init__(self, S0=const.S0, **kwargs):
super(_Insolation, self).__init__(**kwargs)
# initialize diagnostics with correct shape
try:
domain = self.domains['sfc']
except:
domain = self.domains['default']
self.add_diagnostic('insolation', Field(np.zeros(domain.shape), domain=domain))
self.add_diagnostic('coszen', Field(np.zeros(domain.shape), domain=domain))
self.add_diagnostic('irradiance_factor', Field(np.ones(domain.shape), domain=domain))
self.S0 = S0
self.declare_input(['S0'])
@property
def S0(self):
r"""Property of solar constant S0.
The parameter S0 is stored using a python property and can be changed through
``self.S0 = newvalue`` which will also update the parameter dictionary.
.. warning::
changing ``self.param['S0']`` will not work!
:getter: Returns the S0 parameter which is stored in attribute ``self._S0``.
:setter: * sets S0 which is addressed as ``self._S0`` to the new value
* updates the parameter dictionary ``self.param['S0']`` and
* calls method :func:`_compute_fixed`
:type: float
"""
return self._S0
@S0.setter
def S0(self, value):
self._S0 = value
self.param['S0'] = value
self._compute_fixed()
def _coszen_from_insolation(self):
return self.insolation / self.S0
def _compute_fixed(self):
'''Recompute any fixed quantities after a change in parameters'''
pass
def _get_current_insolation(self):
self._compute_fixed()
def _compute(self):
self._get_current_insolation()
return {}
[docs]
class FixedInsolation(_Insolation):
r"""A class for fixed insolation at each point of latitude of the domain.
The solar distribution for the whole domain is constant and specified by
a parameter.
**Initialization parameters**
:param float S0: solar constant
- unit: :math:`\frac{\textrm{W}}{\textrm{m}^2}`
- default value: ``const.S0/4 = 341.2``
:Example:
::
>>> import climlab
>>> from climlab.radiation.insolation import FixedInsolation
>>> model = climlab.EBM()
>>> sfc = model.Ts.domain
>>> fixed_ins = FixedInsolation(S0=340.0, domains=sfc)
>>> print(fixed_ins)
climlab Process of type <class 'climlab.radiation.insolation.FixedInsolation'>.
State variables and domain shapes:
The subprocess tree:
top: <class 'climlab.radiation.insolation.FixedInsolation'>
"""
def __init__(self, S0=const.S0/4, **kwargs):
super(FixedInsolation, self).__init__(S0=S0, **kwargs)
self._compute_fixed()
def _compute_fixed(self):
self.insolation[:] = self.S0
self.coszen[:] = self._coszen_from_insolation()
class _SteadyInsolation(_Insolation):
'''A parent class for processes that use a steady-in-time insolation'''
def _calc_insolation():
# classes need to define the correct method along with necessary input arguments
pass
def _compute_fixed(self):
try:
insolation = self._calc_insolation()
coszen = insolation / self.S0
dom = self.domains['default']
# make sure that the diagnostic has the correct field dimensions.
dom = self.domains['default']
try:
insolation = to_latlon(insolation, domain=dom)
coszen = to_latlon(coszen, domain=dom)
self._steady_insolation = insolation
self._steady_coszen = coszen
self._steady_irradiance_factor = 1.
except:
self._steady_insolation = Field(insolation, domain=dom)
self._steady_coszen = Field(coszen, domain=dom)
self._steady_irradiance_factor = 1.
# Correctly populate the diagnostic fields
self._get_current_insolation()
except AttributeError: # The silent fail is here just for the initialization step
pass
def _get_current_insolation(self):
self.insolation[:] = self._steady_insolation
self.coszen[:] = self._steady_coszen
self.irradiance_factor[:] = self._steady_irradiance_factor
[docs]
class P2Insolation(_SteadyInsolation):
r"""A class for parabolic solar distribution over the domain's latitude
on the basis of the second order Legendre Polynomial.
Calculates the latitude dependent solar distribution as
.. math::
S(\varphi) = \frac{S_0}{4} \left( 1 + s_2 P_2(x) \right)
where :math:`P_2(x) = \frac{1}{2} (3x^2 - 1)` is the second order Legendre Polynomial
and :math:`x=sin(\varphi)`.
**Initialization parameters**
:param float S0: solar constant
- unit: :math:`\frac{\textrm{W}}{\textrm{m}^2}`
- default value: ``1365.2``
:param floar s2: factor for second legendre polynominal term
- default value: ``-0.48``
:Example:
::
>>> import climlab
>>> from climlab.radiation.insolation import P2Insolation
>>> model = climlab.EBM()
>>> sfc = model.Ts.domain
>>> p2_ins = P2Insolation(S0=340.0, s2=-0.5, domains=sfc)
>>> print(p2_ins)
climlab Process of type <class 'climlab.radiation.insolation.P2Insolation'>.
State variables and domain shapes:
The subprocess tree:
top: <class 'climlab.radiation.insolation.P2Insolation'>
"""
def __init__(self, S0=const.S0, s2=-0.48, **kwargs):
super(P2Insolation, self).__init__(S0=S0, **kwargs)
self.s2 = s2
self._compute_fixed()
@property
def s2(self):
r"""Property of second legendre polynomial factor s2.
s2 in following equation:
.. math::
S(\varphi) = \frac{S_0}{4} \left( 1 + s_2 P_2(x) \right)
:getter: Returns the s2 parameter which is stored in attribute ``self._s2``.
:setter: * sets s2 which is addressed as ``self._S0`` to the new value
* updates the parameter dictionary ``self.param['s2']`` and
* calls method :func:`_compute_fixed`
:type: float
"""
return self._s2
@s2.setter
def s2(self, value):
self._s2 = value
self.param['s2'] = value
self._compute_fixed()
def _calc_insolation(self):
phi = np.deg2rad(self.lat)
return self.S0 / 4 * (1. + self.s2 * P2(np.sin(phi)))
# These classes calculate insolation based on orbital parameters
# and astronomical formulas
[docs]
class AnnualMeanInsolation(_SteadyInsolation):
r"""A class for latitudewise solar insolation averaged over a year.
This class computes the daily-mean solar insolation for each day of the year and
latitude specified in the domain on the basis of orbital parameters and
astronomical formulas.
Internally this process calls the method :func:`~climlab.solar.insolation.daily_insolation_factors`
to compute solar zenith angle and adjustments to total irradiance.
See there for details on how the solar distribution depends on orbital parameters.
The mean over the year is calculated from data returned by
:func:`~climlab.solar.insolation.daily_insolation_factors` and stored in the
diagnostics ``insolation``, ``coszen``, and ``irrandiance_factor``.
Different daily averaging methods can be specified for the zenith angle via the
``weighting`` argument. See :func:`~climlab.solar.insolation.daily_insolation_factors`
for more details.
**Initialization parameters** \n
:param float ``S0``: solar constant \n
- unit: :math:`\frac{\textrm{W}}{\textrm{m}^2}` \n
- default value: ``1365.2``
:param dict ``orb``: a dictionary with three orbital parameters (as provided by
:class:`~climlab.solar.orbital.OrbitalTable`):
* ``'ecc'`` - eccentricity
* unit: dimensionless
* default value: ``0.017236``
* ``'long_peri'`` - longitude of perihelion (precession angle)
* unit: degrees
* default value: ``281.37``
* ``'obliquity'`` - obliquity angle
* unit: degrees
* default value: ``23.446``
:param str ``weighting``: flag to specify daily averaging method for solar zenith angle. Valid options are: \n
- ``'time'`` (default): unweighted 24 hour daily average
- ``'sunlit'``: time average over sunlit hours
- ``'insolation'``: insolation-weighted average
**Object attributes** \n
Additional to the parent class :class:`~climlab.radiation.insolation._Insolation`
following object attributes are generated and updated during initialization:
:ivar insolation: Current insolation in W/m2
:vartype insolation: Field
:ivar coszen: Cosine of the current solar zenith angle
:vartype coszen: Field
:ivar dict orb: initialized with given argument ``orb``
:Example:
Create regular EBM and replace standard insolation subprocess by
:class:`~climlab.radation.AnnualMeanInsolation`::
>>> import climlab
>>> from climlab.radiation import AnnualMeanInsolation
>>> # model creation
>>> model = climlab.EBM()
>>> print(model)
.. code-block:: none
:emphasize-lines: 12
climlab Process of type <class 'climlab.model.ebm.EBM'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.AplusBT.AplusBT'>
albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
iceline: <class 'climlab.surface.albedo.Iceline'>
cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.P2Insolation'>
::
>>> # catch model domain for subprocess creation
>>> sfc = model.domains['Ts']
>>> # create AnnualMeanInsolation subprocess
>>> new_insol = AnnualMeanInsolation(domains=sfc, **model.param)
>>> # add it to the model
>>> model.add_subprocess('insolation',new_insol)
>>> print(model)
.. code-block:: none
:emphasize-lines: 12
climlab Process of type <class 'climlab.model.ebm.EBM'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.AplusBT.AplusBT'>
albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
iceline: <class 'climlab.surface.albedo.Iceline'>
cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.AnnualMeanInsolation'>
"""
def __init__(self, S0=const.S0, orb=const.orb_present,
weighting='time', **kwargs):
super(AnnualMeanInsolation, self).__init__(S0=S0, **kwargs)
self.orb = orb
self.weighting = weighting
self._compute_fixed()
@property
def orb(self):
r"""Property of dictionary for orbital parameters.
orb contains: (for more information see :class:`~climlab.solar.orbital.OrbitalTable`)
* ``'ecc'`` - eccentricity [unit: dimensionless]
* ``'long_peri'`` - longitude of perihelion (precession angle) [unit: degrees]
* ``'obliquity'`` - obliquity angle [unit: degrees]
:getter: Returns the orbital dictionary which is stored in attribute
``self._orb``.
:setter: * sets orb which is addressed as ``self._orb`` to the new value
* updates the parameter dictionary ``self.param['orb']`` and
* calls method :func:`_compute_fixed`
:type: dict
"""
return self._orb
@orb.setter
def orb(self, value):
self._orb = value
self.param['orb'] = value
self._compute_fixed()
def _calc_insolation(self):
return annual_insolation(self.lat, orb=self.orb, S0=self.S0, days_per_year=const.days_per_year)
[docs]
class DailyInsolation(AnnualMeanInsolation):
r"""A class to compute latitudewise daily average solar insolation for specific
days of the year.
This class computes the solar insolation on basis of orbital parameters and
astronomical formulas.
Internally this process calls the method :func:`~climlab.solar.insolation.daily_insolation_factors`.
to compute solar zenith angle and adjustments to total irradiance.
See there for details on how the solar distribution depends on orbital parameters.
Different daily averaging methods can be specified for the zenith angle via the
``weighting`` argument. See :func:`~climlab.solar.insolation.daily_insolation_factors`
for more details.
**Initialization parameters**
:param float S0: solar constant
- unit: :math:`\frac{\textrm{W}}{\textrm{m}^2}`
- default value: ``1365.2``
:param dict orb: a dictionary with orbital parameters:
* ``'ecc'`` - eccentricity
* unit: dimensionless
* default value: ``0.017236``
* ``'long_peri'`` - longitude of perihelion (precession angle)
* unit: degrees
* default value: ``281.37``
* ``'obliquity'`` - obliquity angle
* unit: degrees
* default value: ``23.446``
:param str ``weighting``: flag to specify daily averaging method for solar zenith angle. Valid options are: \n
- ``'time'`` (default): unweighted 24 hour daily average
- ``'sunlit'``: time average over sunlit hours
- ``'insolation'``: insolation-weighted average
**Object attributes** \n
Additional to the parent class :class:`~climlab.radiation.insolation._Insolation`
following object attributes are generated and updated during initialization:
:ivar insolation: Current insolation in W/m2
:vartype insolation: Field
:ivar coszen: Cosine of the current solar zenith angle
:vartype coszen: Field
:ivar dict orb: initialized with given argument ``orb``
:Example:
Create regular EBM and replace standard insolation subprocess by
:class:`~climlab.radation.DailyInsolation`::
>>> import climlab
>>> from climlab.radiation import DailyInsolation
>>> # model creation
>>> model = climlab.EBM()
>>> print(model)
.. code-block:: none
:emphasize-lines: 12
climlab Process of type <class 'climlab.model.ebm.EBM'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.AplusBT.AplusBT'>
albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
iceline: <class 'climlab.surface.albedo.Iceline'>
cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.P2Insolation'>
::
>>> # catch model domain for subprocess creation
>>> sfc = model.domains['Ts']
>>> # create DailyInsolation subprocess and add it to the model
>>> model.add_subprocess('insolation',DailyInsolation(domains=sfc, **model.param))
>>> print(model)
.. code-block:: none
:emphasize-lines: 12
climlab Process of type <class 'climlab.model.ebm.EBM'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.AplusBT.AplusBT'>
albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
iceline: <class 'climlab.surface.albedo.Iceline'>
cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.DailyInsolation'>
"""
def _compute_fixed(self):
# One full year of current times
self._times_in_a_year = np.arange(self.current_time,
self.current_time + np.timedelta64(int(const.seconds_per_year), 's'),
self.timestep)
try:
coszen, irradiance_factor = daily_insolation_factors(self.lat,
dates_to_day_index(self._times_in_a_year),
orb=self.orb,
weighting=self.weighting)
dom = self.domains['default']
if 'lon' in dom.axes:
# coszen is latitude-only, need to broadcast across longitude
# assumption is axes are ordered (lat, lon, depth)
# NOTE this is a clunky hack and all this will go away
# when we use xarray structures for these internals
num_lon = dom.axes['lon'].num_points
coszen = np.tile(coszen[...,np.newaxis, :], [1, num_lon, 1])
irradiance_factor = np.tile(irradiance_factor[...,np.newaxis, :], [1, num_lon, 1])
self._coszen_array = coszen
self._irradiance_factor_array = irradiance_factor
except AttributeError:
pass
def _get_current_insolation(self):
now_index = np.where(self._times_in_a_year==self.current_time)[0]
if now_index.size==0:
self._compute_fixed() # No matching date, need to start a new calendar year
now_index = np.where(self._times_in_a_year==self.current_time)[0]
self.coszen[:] = self._coszen_array[..., now_index]
self.irradiance_factor[:] = self._irradiance_factor_array[..., now_index]
self.insolation[:] = self.S0 * self.coszen * self.irradiance_factor
def _compute(self):
self._get_current_insolation()
return {}
[docs]
class InstantInsolation(DailyInsolation):
r"""A class to compute latitudewise instantaneous solar insolation for specific
days of the year.
This class computes the solar insolation on basis of orbital parameters and
astronomical formulas.
Therefore it uses the method :func:`~climlab.solar.insolation.instant_insolation`.
For details how the solar distribution is dependend on orbital parameters
see there.
**Initialization parameters**
:param float S0: solar constant
- unit: :math:`\frac{\textrm{W}}{\textrm{m}^2}`
- default value: ``1365.2``
:param dict orb: a dictionary with orbital parameters:
* ``'ecc'`` - eccentricity
* unit: dimensionless
* default value: ``0.017236``
* ``'long_peri'`` - longitude of perihelion (precession angle)
* unit: degrees
* default value: ``281.37``
* ``'obliquity'`` - obliquity angle
* unit: degrees
* default value: ``23.446``
**Object attributes** \n
Additional to the parent class :class:`~climlab.radiation.insolation._Insolation`
following object attributes are generated and updated during initialization:
:ivar insolation: Current insolation in W/m2
:vartype insolation: Field
:ivar coszen: Cosine of the current solar zenith angle
:vartype coszen: Field
:ivar dict orb: initialized with given argument ``orb``
:Example:
Create regular EBM and replace standard insolation subprocess by
:class:`~climlab.radation.DailyInsolation`::
>>> import climlab
>>> from climlab.radiation import InstantInsolation
>>> # model creation
>>> model = climlab.EBM()
>>> print(model)
.. code-block:: none
:emphasize-lines: 12
climlab Process of type <class 'climlab.model.ebm.EBM'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.AplusBT.AplusBT'>
albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
iceline: <class 'climlab.surface.albedo.Iceline'>
cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.P2Insolation'>
::
>>> # catch model domain for subprocess creation
>>> sfc = model.domains['Ts']
>>> # create InstantInsolation subprocess and add it to the model
>>> model.add_subprocess('insolation',InstantInsolation(domains=sfc, **model.param))
>>> print(model)
.. code-block:: none
:emphasize-lines: 12
climlab Process of type <class 'climlab.model.ebm.EBM'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.AplusBT.AplusBT'>
albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
iceline: <class 'climlab.surface.albedo.Iceline'>
cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.InstantInsolation'>
"""
def _compute_fixed(self):
# One full year of current times
self._times_in_a_year = np.arange(self.current_time,
self.current_time + np.timedelta64(int(const.seconds_per_year), 's'),
self.timestep)
try:
dom = self.domains['default']
if 'lon' in dom.axes:
lon = self.lon
else:
lon = 0.
coszen, irradiance_factor = instant_insolation_factors(self.lat,
dates_to_day_index(self._times_in_a_year),
lon=lon, orb=self.orb)
self._coszen_array = coszen
self._irradiance_factor_array = irradiance_factor
except AttributeError:
pass
def _get_current_insolation(self):
# make sure that the diagnostic has the correct field dimensions.
dom = self.domains['default']
lon = 0
if 'lon' in dom.axes:
lon = self.lon
coszen, irradiance_factor = instant_insolation_factors(self.lat,
dates_to_day_index(self.current_time),
lon=lon, orb=self.orb,)
self.coszen[:] = Field(coszen, domain=dom)
self.irradiance_factor[:] = Field(irradiance_factor, domain=dom)
self.insolation[:] = self.S0 * self.coszen * self.irradiance_factor