daily_insolation documentation

daily_insolation computes daily average insolation as a function of day and latitude at any point during the past 5 million years.

This function is from Ian Eisenman and Peter Huybers (see the reference below). The daily_insolation function is quite similar to the solar_radiation function, but one may suit your needs better than the other. The daily_insolation function is best suited for investigations involving orbital changes over thousands to millions of years, whereas the solar_radiation may be easier to use for applications such as present-day precipitation/drought research.

Back to Climate Data Tools Contents.

Contents

Syntax

Fsw = daily_insolation(kyear,lat,day)
Fsw = daily_insolation(kyear,lat,day,day_type)
Fsw = daily_insolation(kyear,lat,day,day_type,'constant',So)
Fsw = daily_insolation(kyear,lat,day,day_type,'mjmd')
[Fsw, ecc, obliquity, long_perh] = daily_insolation(...)

Description

Fsw = daily_insolation(kyear,lat,day) gives the daily average insolation in W/m^2 at latitude(s) on specified day(s) of the year (as given by the doy function) for the times kyear, which are thousands of years before present. For example, use kyear = +3000 to indicate 3 million years before present. Maximum allowed value of kyear is 5000.

Fsw = daily_insolation(kyear,lat,day,day_type) specifies an optional day type as 1 (default) or 2. The default option 1 specifies days in the range 1 to 365.24, where day 1 is January first and the vernal equinox always occurs on day 80. Option 2 specifies day input as solar longitude in the range 0 to 360 degrees. Solar longitude is the angle of the Earth's orbit measured from spring equinox (21 March). Note that calendar days and solar longitude are not linearly related because, by Kepler's Second Law, Earth's angular velocity varies according to its distance from the sun. If day_type is negative, kyear is taken to be a 3 element array containing [eccentricity, obliquity, and longitude of perihelion].

Fsw = daily_insolation(kyear,lat,day,day_type,'constant',So) specfies a solar constant So. Default So is 1365 W/m^2.

Fsw = daily_insolation(kyear,lat,day,day_type,'mjmd') returns Fsw in units of (MJ/m^2)/day rather than the default W/m^2.

[Fsw, ecc, obliquity, long_perh] = daily_insolation(...) also returns the orbital eccentricity, obliquity, and longitude of perihelion (precession angle).

Example 1

For this example, calculate summer solstice insolation at 65 N. We'll do it two different ways. First, use doy to get the day of year corresponding to solstice, assuming it occurs on June 20:

kyr = 0:1000;

Fsw_jun20 = daily_insolation(kyr,65,doy('june 20'));

plot(kyr,Fsw_jun20)
set(gca,'xdir','reverse') % reverses x direction
ylabel('summer solstice insolation at 65\circN (W m^{-1})')
xlabel 'thousands of years before present'

Rather than approximating June 20 as the solstice, a more precise way would be to specify the exact orbital angle corresponding to the solstice (90 degrees) using the day_type=2 option, like this:

Fsw_exact = daily_insolation(kyr,65,90,2);

hold on
plot(kyr,Fsw_exact)

legend('June 20','exact')

Example 2

We can plot the difference between June 20 (calendar day) and the exact summer solstice insolation at 65 N

Fsw_jun20 = daily_insolation(kyr,65,doy('june 20'));
Fsw_exact = daily_insolation(kyr,65,90,2);

figure
plot(kyr,Fsw_jun20-Fsw_exact)
set(gca,'xdir','reverse') % reverses x direction
ylabel('difference in solstice insolation at 65\circN (W m^{-1})')
xlabel 'thousands of years before present'

Example 3

Insolation for the current orbital configuration as a function of day and latitude. Since this is a function of two variables, use meshgrid to create grids of the day and latitude variables:

[day,lat]=meshgrid(1:5:365, -90:90);
[Fsw,ecc,obl,omega]=daily_insolation(0,lat,day);

disp([ecc,obl,omega]) % displays optional outputs
    0.0172   23.4460  281.3700

We can now plot insolation as a function of day and latitude like this:

figure
[c,h]=contour(day,lat,Fsw,[0:50:500]);
clabel(c,h)
xlabel 'day of year'
ylabel 'latitude'
title 'present day insolation (W m^{-2})'

Example 4

Calculate annual (not daily) average insolation by explicitly specifying orbital parameters.

ecc=0.017236;
obl=23.446;
omega=101.37+180;
[day,lat]=meshgrid(1:5:365, -90:90);

Fsw=daily_insolation([ecc,obl,omega],lat,day,-1);

Fsw_annual = mean(Fsw,2);


figure
plot(-90:90,Fsw_annual)
axis tight
xlabel 'latitude'
ylabel 'annual mean insolation (W m^{-2})'

Example 5

Compare calculated insolation with example values given by Berger (1991). Start by loading Earth's orbital parameters from Berger A. and Loutre M.F., 1991, which is included as a sample dataset in CDT. Note, Berger and Loutre used negative values for their kiloyears, so we'll have to multiply them by -1. Also, they used a solar constant of 1360 W/m^2 so we'll specify that as well.

D = importdata('orbit91.txt'); %

kyear = -D.data(:,1);
insol_65NJul = D.data(:,6);

Fsw=daily_insolation(kyear,65,(7-3)*30,2,'constant',1360);

figure
plot(kyear,1-insol_65NJul./Fsw)

Note, the values above agree within 3e-5.

Example 6

Plot 65N integrated summer insolation as in Huybers (2006), Science 313 508-511:

[kyear,day]=meshgrid(1000:1:2000,1:1:365);
Fsw = daily_insolation(kyear,65,day);
Fsw(Fsw<275)=0;

figure
plot(-(1000:1:2000),sum(Fsw,1)*86400*10^-9)
title('As in Huybers (2006) Fig. 2C')

Detailed description of calculation:

Values for eccentricity, obliquity, and longitude of perihelion for the past 5 Myr are taken from Berger and Loutre 1991 (data from ncdc.noaa.gov). If using calendar days, solar longitude is found using an approximate solution to the differential equation representing conservation of angular momentum (Kepler's Second Law). Given the orbital parameters and solar longitude, daily average insolation is calculated exactly following Berger 1978.

References:

Authors:

Ian Eisenman and Peter Huybers, Harvard University, August 2006 eisenman@fas.harvard.edu. Slightly modified by Chad A. Greene for inclusion in the Climate Data Toolbox for Matlab, 2019.