VI/63 Earth Orbit, Precession and Insolation -20Myr to +10Myr (Laskar+ 1993)
Orbital, precessional, and insolation quantities for the Earth from
-20 Myr to +10Myr
Laskar J., Joutel F., Boudin F.
<Astron. Astrophys. 270, 522 (1993)>
=1993A&A...270..522L 1993A&A...270..522L
ADC_Keywords: Earth ; Solar system
Keywords: celestial mechanics: insolation, precession - Earth -
solar system: numerical methods, resonances
********************************************************************************
* *
* NOTICE: *
* *
* This software can be freely copied and distributed provided that: *
* * All the included files are kept together. *
* * This notice is included in the distribution *
* * no changes or alterations are made on the distributed files *
* *
* Any changes can be performed on the FORTRAN files, for specific use, *
* but these changes cannot be included in the distribution files, and *
* reference to the original programs should be included. *
* *
* WARNING: *
* This software is not designed to be used in any possible conditions, *
* and is not presently in a definite version. *
* *
* In particular, the user must take care that the starting and ending *
* times for which he computes data are contained in the necessary files. *
* The stepsize should also be the same. *
* *
* THE OBJECT OF THIS WORK IS TO PROVIDE A USEFUL TOOL FOR THE PALEOCLIMATE *
* COMMUNITY. ANY COMMENTS OR WISH FOR IMPROVEMENTS ARE WELCOME. *
* *
* Authors: J. Laskar, F. Joutel and some other contributions *
* (c) Astronomie et Systemes Dynamiques, Bureau des Longitudes, Paris (1993) *
* *
* Jacques Laskar *
* Astronomie et Systemes Dynamiques, *
* Bureau des Longitudes *
* 77 av. Denfert-Rochereau *
* 75014 Paris *
* email: *
* *
********************************************************************************
File Summary:
(the second part lists files which can be generated in the nominal solution)
--------------------------------------------------------------------------------
FileName Lrecl Records Explanations
--------------------------------------------------------------------------------
ReadMe 80 . This file
Makefile 80 98 to compile and generate derived files
ORBELN.ASC 107 20101 Normal orbital solution La93(0,1),
before present years (-20Myr to 0)
ORBELP.ASC 107 10101 Normal orbital solution La93(0,1)
after present years (0 to +10Myr)
climavar.f 75 86 FORTRAN program for computation of usual
climatic variables
climavar.par 28 11 parameter file for climavar
insola.f 79 253 FORTRAN program for computations of various
insolation quantities
insola.par 24 8 parameter file for insola
insolsub.f 79 433 subroutines for insola.f
integ.f 75 300 FORTRAN program for integration step
integ.par 26 13 parameter file for integ
integsub.f 82 882 FORTRAN subroutine for integ.f
prepa.f 80 709 FORTRAN program for preparation step
prepa.par 26 9 parameter file for prepa
prepinsol.f 75 168 FORTRAN program for preparation step for insol
prepinsol.par 26 10 parameter file for prepinsol
CLIVAR0N.ASC 71 20001 Usual climatic variables in Nominal solution
La93(0,1) before present years (-20Myr to 0)
CLIVAR0P.ASC 71 10001 Usual climatic variables in Nominal solution
La93(0,1) after present years (0 to +10Myr)
PREC0N.ASC 59 20051 Nominal solution La93(0,1)
before present years (-20Myr to 0)
PREC0P.ASC 59 10051 Nominal solution La93(0,1)
after present years (0 to +10Myr)
--------------------------------------------------------------------------------
Byte-per-byte Description of file: ORBELN.ASC ORBELP.ASC
--------------------------------------------------------------------------------
Bytes Format Units Label Explanations
--------------------------------------------------------------------------------
4- 11 F8.0 1000yr t Time from J2000 in 1000 years
13- 35 D23.16 --- k = e * cos (longitude of perihelion)
(e: eccentricity)
37- 59 D23.16 --- h = e * sin (longitude of perihelion)
61- 83 D23.16 --- q = sin(i/2) cos (longitude of node)
(i: inclination from J2000 ecliptic)
85-107 D23.16 --- p = sin(i/2) sin (longitude of node)
--------------------------------------------------------------------------------
Byte-per-byte Description of file: PREC0N.ASC PREC0P.ASC
--------------------------------------------------------------------------------
Bytes Format Units Label Explanations
--------------------------------------------------------------------------------
4- 11 F8.0 1000yr t Time from J2000 in 1000 years
13- 35 D23.16 rad eps obliquity (radians)
37- 59 D23.16 rad phi general precession in longitude (radians)
--------------------------------------------------------------------------------
Byte-per-byte Description of file: CLIVAR0N.ASC CLIVAR0P.ASC
--------------------------------------------------------------------------------
Bytes Format Units Label Explanations
--------------------------------------------------------------------------------
4- 11 F8.0 1000yr t Time from J2000 in 1000 years
13- 31 D19.12 --- e Eccentricity
33- 51 D19.12 rad eps obliquity (radians)
53- 71 D19.12 --- CP e * sin(longitude of perihelion
from moving equinox)
--------------------------------------------------------------------------------
Description:
La93 is a program for computing the precession and obliquity of
the Earth for various values of
1) the tidal effect of the Moon (CMAR)
2) the dynamical ellipticity of the Earth (FGAM)
The nominal solution La90 corresponds to (CMAR = 0., FGAM = 1).
The general solution will be called La93(CMAR,FGAM), thus
La90 = La93(0.,1.)
and La93(1.,1.) is obtained with the same tidal effect as
in Quinn, Tremaine, Duncan (1991), although these solutions
are not completely identical (see Laskar, Joutel, Boudin, 1993)
The files and software of this package can be used in three
different manners:
1) Contruction and Use of the nominal solution La93(0,1)
The ASCII files ORBEL*.ASC contain the nominal orbital solution.
The ASCII files PREC0*.ASC contain the nominal precession solution.
The ASCII files CLIVAR0*.ASC contain the nominal climatic solution.
The files PREC*.ASC and CLIVAR*.ASC can also be generated from
the enclosed files (see section 2)
For the computation of insolation quantities, the user will execute
the 'prepinsol' step, and then 'insola'.
2) Construction of a parametrized La93(CMAR, FGAM) new solution
The user reconstructs a complete La93(CMAR, FGAM) solution.
The compilation of all required programs is obtained by running
the command 'make' on a Unix machine.
The preparation step 'prepa' needs to be done once, in order
to prepare the necessary binary files.
Then 'integ' will construct the new solutions for the given
parameters (CMAR, FGAM).
Alternatively, change in the Makefile the values of CMAR and FGAM
before running the 'make clean' command (removes the files computed
using the preceding values of CMAR and FGAM), and 'make La93'.
3) Changes in the model of precession, for example to take into account
some feedback resulting from redistribution of the ice on the Earth
resulting from climate changes.
In this case, and in this case only, the user needs to edit
the FORTRAN file integ.f
More precisely, the subroutines which can be eventually
modified are
SUBROUTINE INIPRE(IPT)
SUBROUTINE PRECES(t,AK,AH,AQ,AP,DK,DH,DQ,DP,AKI,DKI)
The users may also want to adapt the driver INTEG to his specific needs
References:
La90: Laskar, J.: 1990, The chaotic motion of the solar system.
A numerical estimate of the chaotic zones
Icarus, 88, 266
La93: Laskar, J., Joutel, F., Boudin, F.: 1993, Orbital, precessional
and insolation quantities for the Earth
from -20 Myr to + 10Myr
Astron. Astrophys. 270, 522
Detailed Description of the Program
-----------------------------------
The program consists of several steps
prepa : preparation step.
This step generates the necessary binary files
from the ASCII files containing the elliptical elements
of the Earth every 1000 yr, from -20Myr to +10 Myr
from the solution La90.
THE USER should edit the file of command prepa.par
run prepa
two new files : nomfich and nomficder will be created
containing the necessary data for the remaining steps
for the time span from datedebut to datefin.
(a user can choose to work only from -10Myr to 0)
integ : integration step.
This step is the integration of the equations of precession
for various values of the tidal effect, and of the dynamical
ellipticity of the Earth.
THE USER should edit the file of command integ.par
run integ
THE USER will be prompted for the values of CMAR and FGAM
during the execution of the program
climavar : computation of the usual climatic variables
prepinsol : preparation of the insolation computations
insola : interactive insolation computations
********************************************************************
* NOTE: in prepa.par, integ.par, and prepinsol.par, *
* the starting date 'datedebut' and final date 'datefin' *
* need to be compatible, otherwise unexpected behaviour will occur.*
* In the distributed files, these are set * to -20 Myr and 10 Myr, *
* in order to obtain the examples ASCII files of the *
* nominal solution. But the user can reduce this to a much shorter *
* interval in order to gain computing time and disk space. *
********************************************************************
--------------------------------------------------------------------------------
PREPA
This step uses the parameter file prepa.par with the following meaning:
nomascpos: ASCII file for positive time (IN)
nomascneg: ASCII file for negative time (IN)
nomfich : Binary file for elliptical elements (OUT)
nomficder: Binary file for derivatives (OUT)
datedebut: starting time (Myr)
datefin : ending time (Myr)
( -20 ≤ datedebut < datefin ≤ +10 )
statut : 'new' or 'unknown' status of the created files
Example of prepa.par:
--------------------------
&NAMSTD
nomfich = 'ELL.BIN',
nomfichder = 'DER.BIN',
nomascpos = 'ORBELP.ASC',
nomascneg = 'ORBELN.ASC',
datedebut=-20.D0,
datefin=10.D0,
statut='unknown'
&END
--------------------------
The user may restrict himself to e.g. -2Myr to 0 by putting instead
datedebut=-2.D0
datefin=0.D0
The output of the 'prepa' command gives the following ouput:
-------------------------------------------------------------------------
Solution La90-La93 for the Earth
La93.prepa version 0.8
preparation step
(c) ASD/BdL (1993)
ASCII file for positive time :
ORBELP.ASC
ASCII file for negative time :
ORBELN.ASC
Binary file for elliptical elements :
ELL.BIN
Binary file for derivatives of ell. el :
DER.BIN
starting time (Myr) : -20.00000000000000
ending time (Myr) : 0.0000000000000000E+00
The file ELL.BIN is created.
The file DER.BIN is created.
The preparation step is completed normally
-------------------------------------------------------------------------
(Duration: about 15 minutes are necessary on a Sparc2 station)
--------------------------------------------------------------------------------
INTEG
This step uses the parameter file integ.par with the following meaning:
nomfich : binary file for the elements t,k,h,q,p
this file is computed in the preparation phase
prepa. (IN)
nomfichder : binary file for the derivatives of the elliptical
elements. This file is computed in the preparation
phase prepa. (IN)
pas : integration step (in years). The default, and
advised value is 200
nechant : sampling step. The results are written in
the files every nechant*pas years.
default nechant = 5
the files gives positions every 5*200 = 1000 years
datefin : ending time (Myr) (≥0)
datedebut : starting time (Myr) (≤0)
ecritpos : 'oui' for output in an ASCII file for positive times
'non' screen output
ecritneg : 'oui' for output in an ASCII file for negative times
'non' screen output
statut : status of the created files ( 'new' or 'unknown' )
fichrespos : name of ASCII file for positif times
t,eps,phi (OUT)
t : time (unit = 1000yr)
eps : obliquity (radians)
phi : general precession in longitude (radians)
fichresneg : name of ASCII file for negative times
t,eps,phi (OUT)
t : time (unit = 1000yr)
eps : obliquity (radians)
phi : general precession in longitude (radians)
Example of integ.par
--------------------------
&NAMSTD
nomfich = 'ELL.BIN',
nomfichder = 'DER.BIN',
pas = 200,
nechant=5,
datefin= 10.D0,
datedebut= -20.D0,
statut= 'unknown',
ecritpos = 'oui',
ecritneg = 'oui',
fichrespos = 'PRECP.ASC',
fichresneg = 'PRECN.ASC'
&END
--------------------------
The dialog of the 'integ' command gives the following ouput:
-------------------------------------------------------------------------
La93.integ version 0.8
integration step
(c) ASD/BdL (1993)
stepsize (yrs) default : 200 : 200.0000000000000
print every 5 step
starting time (Myr) : -5.000000000000000
ending time (Myr) : 0.0000000000000000E+00
save results for positive time :oui
save results for negative time :oui
ASCII file for positive time :
PRECP.ASC
ASCII file for negative time :
PRECN.ASC
status of created files :unknown
Relative change of dynamical ellipticity
ENTER gamma/gamma_0 (default 1)
1 ≪≪≪----------------------------------- the user need to answer
FGAM = 1.000000000000000
Relative tidal effect
0 : no tidal effect
1 : -4.6 D-18 seconds**-1
ENTER cmar (default 0)
0 ≪≪≪----------------------------------- the user need to answer
CMAR = 0.0000000000000000E+00
Informations for internal check
VITESSE ANGULAIRE DE LA TERRE 474659981.5971373
PRECESSION EN ASCENSION DROITE 5038.783333037391
ELD 3.2800511416867092E-03
RFL0 34.42998531044029
RFL1 -2.6896187697037889E-03
RFL3 3.3459472599329403E-04
RFS 15.97940304351487
RFL0+RFL1+RFL3+RFS 50.40703332991146
CP1,CP2,CP3,CP4 EN SECONDES PAR AN
37.52660322621580 -1.5651731683509029E-03 8.2602928316139284E-05
34.81861759592061
AK1 0.0000000000000000E+00
AK2 0.0000000000000000E+00
-------------------------------------------------------------------------
(Duration: about 15 minutes are necessary on a Sparc2 station)
************************************************************
* REMARK: when one ask for negative and positive time, *
* the user is prompted twice for the preceding questions *
************************************************************
--------------------------------------------------------------------------------
CLIMAVAR
After this integration step, if the user prefers the usual
climatic variables, he can run the small program
climavar
which is provided as a FORTRAN source code;
the parameters of climavar are taken in climavar.par with similar
definitions as in the preceding steps.
Alternatively, one can execute
make climate
--------------------------------------------------------------------------------
INSOLA
The two routines prepinsol and insola are designed to compute all
necessary insolation quantities derived from the orbital and
precessional quantities computed above.
They are given on the form of FORTRAN code, so the user can check
if they correspond to his needs. He can also design his own insolation
routines.
A preparation step prepinsol is necessary prior to execution of
'insola'. This preparation step will construct a binary file with the
useful quantities.
a) set the parameters of prepinsol.par
Example:
---------------------------
&NAMSTD
nomascpos = 'ORBELP.ASC',
nomascneg = 'ORBELN.ASC',
nomprecpos = 'PRECP.ASC',
nomprecneg = 'PRECN.ASC',
nomfich = 'SOLCLI.BIN',
datedebut= -20.D0,
datefin= 10.D0,
statut='unknown'
&END
---------------------------
b) run prepinsol
(these steps are automatically done with the 'make climate' command)
For the actual computations of the insolation parameters, the user has to:
a) set the parameters of insola.par; in this file, the 'so' parameter
refers to the Solar constant expressed in W/m2.
Example:
------------------------
&NAMSTD
nomfich = 'SOLCLI.BIN',
pas = 1.D3,
datedebut = -1.D0,
datefin =0.D0,
statut='unknown',
so = 1350.D0,
&END
------------------------
b) run insola
insola is a self documented program;
For more details, the user can refer to La93 (reference above)
(End) Jacques Laskar [Bureau des Longitudes] 05-Nov-1994