Insol_l1l2: Time-integrated insolation

In mcrucifix/palinsol: Insolation for Palaeoclimate Studies

Time-integrated insolation

Description

Computes time-integrated incoming solar radiation (Insol) either between given true solar longitudes (Insol_l1l2) or days of year (Insol_d1d2) for a given orbit and latitude

Usage

Insol_l1l2(
  orbit,
  l1 = 0,
  l2 = 2 * pi,
  lat = 65 * pi/180,
  avg = FALSE,
  ell = TRUE,
  ...
)

Insol_d1d2(orbit, d1, d2, lat = 65 * pi/180, avg = FALSE, ...)

Arguments

orbit	Output from a solution, such as ber78, ber90 or la04
l1	lower true solar longitude bound of the time-integral
l2	upper true solar longitude bound of the time-integral
lat	latitude
avg	performs a time-average.
ell	uses elliptic integrals for the calculation (much faster)
...	other arguments to be passed to Insol
d1	lower calendar day (360-day year) of the time-integral
d2	upper calendar day (360-day year) of the time-integral

Details

All angles input measured in radians.

Note that in contrast to Berger (2010) we consider the tropic year as the reference, rather than the sideral year, which partly explains some of the small differences with the original publication

Value

Time-integrated insolation in kJ/m2 if avg=TRUE, else time-average in W/m2

Functions

• Insol_d1d2(): Mean and integrated insolation over an interval bounded by calendar days

Author(s)

Michel Crucifix, U. catholique de Louvain, Belgium.

References

Berger, A., Loutre, M.F. and Yin Q. (2010), Total irradiation during any time interval of the year using elliptic integrals, Quaternary Science Reviews, 29, 1968 - 1982, doi:10.1016/j.quascirev.2010.05.007

Examples

## reproduces Table 1a of Berger et al. 2010:
lat <- seq(85, -85, -10) * pi/180. ## angles in radians. 
orbit=c(eps=  23.446 * pi/180., ecc= 0.016724, varpi= (102.04 - 180)* pi/180. )
T <-  sapply(lat, function(x) c(lat = x * 180/pi, 
        m1 =  Insol_l1l2(orbit, 0, 70 * pi/180, lat=x, ell= TRUE, S0=1368) / 1e3,
        m2 =  Insol_l1l2(orbit, 0, 70 * pi/180, lat=x, ell=FALSE, S0=1368) / 1e3) ) 
data.frame(t(T))
 ## reproduces Table 1b of Berger et al. 2010:
lat <- c(85, 55, 0, -55, -85) * pi/180. ## angles in radians. 
T <-  sapply(lat, function(x) c(lat = x * 180/pi, 
         m1 =  Insol_l1l2(orbit, 30 * pi/180. , 75 * pi/180, 
               lat=x, ell= TRUE, S0=1368) / 1e3,
         m2 =  Insol_l1l2(orbit, 30 * pi/180. , 75 * pi/180, 
               lat=x, ell=FALSE, S0=1368) / 1e3) ) 
 ## reproduces Table 2a of Berger et al. 2010:
lat <- seq(85, -85, -10) * pi/180. ## angles in radians. 

## 21 march in a 360-d year. By definition : day 80 = 21 march at 12u
d1 = 79.5 
d2 = 79.5 + (10 + 30 + 30 ) * 360/365.2425 ## 30th May in a 360-d year

T <-  sapply(lat, function(x) c(lat = x * 180/pi, 
        m1 =  Insol_d1d2(orbit, d1,d2, lat=x, ell= TRUE, S0=1368) / 1e3,
        m2 =  Insol_d1d2(orbit, d1,d2, lat=x, ell= FALSE, S0=1368) / 1e3))
                          
## I did not quite get the same results as on the table 
## on this one; probably a matter of calendar
## note : the authors in fact used S0=1368 (pers. comm.) 
## 1366 in the paper is a misprint

data.frame(t(T))

## annual mean insolation at 65N North, as a function of the longitude of the perihelion
## (expected to be invariant)

varpis <- seq(0,360,10)*pi/180.  
sapply(varpis, function(varpi)
   {   orbit=c(eps=  23.446 * pi/180., ecc= 0.016724, varpi= varpi )
       amean <- Insol_l1l2 (orbit, lat=65*pi/180., avg=TRUE)
       return(amean)
   })

mcrucifix/palinsol documentation built on Oct. 24, 2023, 1:35 a.m.