program integ
C**************************************************************************
c
c                     Driver of the integration
c                      
c
c  Integration of the precession equations are made with the
c   Adams  predictor-corrector method.
c 
c  the following parameters are taken in the  NAMELIST :
c  'integration.par'
c
c           nomfich     : fichier binaire des elements t,k,h,q,p
c           nomfichder  : fichier binaire des derivees des elements
c           pas         : le pas d'integration en annees
c           nechant     : la frequence d'inscription, les resultats
c                          etant inscrits tous les nechant*pas annees
c                          REMARQUE: LE PAS D'ECHANTILLONAGE POUR LA
c                          CREATION DES FICHIERS DESTINES A L'INSOLATION 
c                          EST DONC EGAL A NECHANT*PAS ANNEES         
c           datefin     : la date finale en millions d'annees (>=0)
c           datedebut   : la date initiale en millions d'annees (<=0)
c           ecritpos    : 'oui' si l'on desire l'ecriture des resultats
c                         pour les temps positifs dans un fichier ASCII,
c                         'non' sinon (sortie ‚cran)
c           ecritneg    : 'oui' si l'on desire l'ecriture des resultats
c                         pour les temps negatifs dans un fichier ASCII,
c                         'non' sinon (sortie ‚cran)
c           statut      : pour l'ouverture des nouveaux fichiers
c                         'new' ou 'unknown' 
c           fichrespos  : fichier ASCII pour les resultats t,eps,phi
c                          pour les temps positifs. 
c           fichresneg  : fichier ASCII pour les resultats t,eps,phi
c                          pour les temps n‚gatifs.
c
c   Les autres parametres:
c
c           nbf         : 2 x 4
c           nbel1       : 4 + 1
c           nacd        : nombre de n-uplets dans un enregistrement
c
c  (c) Astronomie et Systemes Dynamiques/Bureau des Longitudes (1993)
c  J. Laskar, F. Joutel 
c  version.0.8 (6/6/1993)
c**************************************************************************
c
      implicit double precision (a-h,o-z)
      parameter (nbf=8,nbel1=5,nacd=100)
      character*3 ecriture,ecritpos,ecritneg
      character*20 statut
      character*50 nomfich,nomfichder,fichrespos,fichresneg
      common/fichier/nechant,ecriture
      common/date/itdebut
      NAMELIST/NAMSTD/nomfich,nomfichder,pas,datefin,datedebut,
     &ecritpos,ecritneg,fichrespos,fichresneg,statut,
     &nechant
      open(8,file='integ.par',form='formatted',status='old')
      read(8,NAMSTD)
      close(8)
      write(*,*) ' Solution La90-La93 for the Earth '
      write(*,*) 
      write(*,*) ' La93.integ  version 0.82'
      write(*,*) ' integration step '
      write(*,*) ' (c) ASD/BdL (1993) '
      write(*,*) 
      write(*,*) ' stepsize (yrs) default : 200           :', pas
      write(*,*) ' print every ',nechant,'  step'
      write(*,*) ' starting time (Myr)                    :',datedebut
      write(*,*) ' ending   time (Myr)                    :',datefin
      write(*,*)
      write(*,*) ' save results for positive time         :',ecritpos
      write(*,*) ' save results for negative time         :',ecritneg
      if (ecritpos.eq.'oui') then
      write(*,*) ' ASCII file for positive time           :',fichrespos
      endif
      if (ecritneg.eq.'oui') then
      write(*,*) ' ASCII file for negative time           :',fichresneg
      endif
      write(*,*) ' status of created files                :', statut
      write(*,*)
      write(*,*)
c
c Les fichiers ont une longueur d'enregistrement correspondant
c a 100 000 ans.
c
      nrecl=8*nbel1*nacd
      open(10,file=nomfich,status='old',access='direct',recl=nrecl)
      open(11,file=nomfichder,status='old',access='direct',recl=nrecl)
      read(10,rec=2) debut
c on possede un un enregistrement de plus aux bornes
      debut=debut-nacd
      itdebut=int(debut)
      if (datefin.GT.0.D0) then
        if (ecritpos.eq.'oui') then
           open(98,file=fichrespos,status=statut)
        endif
        ecriture=ecritpos
        call integration(datefin,pas)
        if (ecriture.eq.'oui') then
            close(98)
        endif
      endif
      if (datedebut .lt. 0D0) then
        if (ecritneg.eq.'oui') then
            open(98,file=fichresneg,status=statut)
        endif
        ecriture=ecritneg
        call integration(datedebut,pas)
        if (ecriture.eq.'oui') then
           close(98)
        endif
      endif
      end

      subroutine integration(t1,pas)
c----------------------------------------------------------------------
c   integration entre 0 et t1 (en millions d'ann‚es) 
c   avec un pas pas (en ann‚es)
c  
c   Routines externes
c           sortie des r‚sultats       :  solout (connection 98)
c           calcul des seconds membres : fcn
c 
c   Parametres:
c
c           nd          : la dimension du systŠme diff‚rentiel (nd = 2)
c           m           : l'ordre de la m‚thode d'Adams
c           ni          : le nombre de points n‚cessaires …
c                         l'interpolation des donn‚es
c
c----------------------------------------------------------------------

      implicit double precision (a-h,o-z)
      external fcn,solout
      parameter (nd=2,m=12,ni=8)
      parameter (nacd=100)
      dimension aki(nd),valyp(nd,0:m-1)
      common/cpres0/akix0,akiy0
      common/indice/ind
      ind=0
c  t1 en millions d'ann‚es
      pas=sign(pas,t1)
c Afin d'avoir suffisemment de donnees pour la suite:
      tfin=t1*1.d6+sign(1.D0,pas)*(nacd/2)*1.d3
c
c Initiatialisation pour la precession:
      IPT=1
      call INIPRE(IPT)
      t0=0.0D0
      aki(1)=akix0
      aki(2)=akiy0

c
c------------------------------------------------------------------------
c Depart de la methode de Runge-Kutta (7)8
c------------------------------------------------------------------------
c
      call depart(nd,fcn,t0,pas,m,aki,valyp,solout)
c
c------------------------------------------------------------------------
c Integration (Adams)
c------------------------------------------------------------------------
c
      call adams(nd,fcn,t0,tfin,pas,aki,valyp,solout)
      return
      end

      SUBROUTINE INIPRE(IPT)
***********************************************************************
*                                                                     *
*     INIPRE  INITIALISATION DE LA PRECESSION POUR LA TERRE           *
*  (c) Astronomie et Systemes Dynamiques, Bureau des Longitudes (1993)*
***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Y)
      DIMENSION aki(2),Dki(2)
      LOGICAL MAREE
      DATA MAREE/.TRUE./
      COMMON/CPRES0/AKIX0,AKIY0
      COMMON/CONST/APG,CP1,CP2,CP3,CP4,AK1,AK2
      AKIY0=0.D0
      RM0=496303.3D-6
      RM1=-20.7D-6
      RM2=-0.1D-6
      RM3=3020.2D-6
      RS0=500209.034508230784D-6
      RS2=0.2D-6
      RDSC=206264.80624709D0
      EPS0=23*3600+26*60+21.448D0
      EPS0=EPS0/RDSC
      AKIX0=DSIN(EPS0)
      GK=0.01720209895D0*365.25D0
      TAML=1/270 68736.47 D0
      TAMT=1/332 946.0D0
      TAMS=1.D0
      AL=384 747.981D0/1.495 978 701D8
      AS=1.000 001 017 78 D0
      OM=72.921151467D-6*86400*365.25
      AMN=-190.771235D0*365.25/RDSC
      AML=17325593.437360D0/RDSC
      BP0=101 803 910 D-15
      APUAI=5029.0966D0/RDSC/100
      APP=APUAI+2*BP0/DTAN(EPS0)
      APG=1.92D0/RDSC/100
      APLS=APP+APG
      RKL=3*GK**2*TAML/AL**3/OM
      RKS=3*GK**2*TAMS/AS**3/OM
      XL0=(RM0-RM2/2)*DCOS(EPS0)
      XL1=RM1*DCOS(2*EPS0)/DSIN(EPS0)
      XL3=-RM3*TAML/(TAML+TAMT)*AML**2/AMN/OM*(6*DCOS(EPS0)**2-1)
      XS=(RS0-RS2/2)*DCOS(EPS0)
      AA=RKL*XL3
      BB=RKL*(XL0+XL1)+RKS*XS
      CC=-APLS
      ELD=-BB+DSQRT(BB**2-4*AA*CC)
      ELD=ELD/2/AA
C88888888
c      ELD=0.32739935D-2
      RKS=3*GK**2*TAMS/AS**3/OM
      RFL0=RKL*ELD*XL0*RDSC
      RFL1=RKL*ELD*XL1*RDSC
      RFL3=RKL*ELD**2*XL3*RDSC
      RFS=RKS*ELD*XS*RDSC
      CP1=RKL*ELD*(RM0-RM2/2)
      CP2=RKL*ELD*RM1
      CP3=-RKL*ELD**2*RM3*TAML/(TAML+TAMT)*AML**2/AMN/OM
      CP4=RKS*ELD
***-------------------- Changement d'ellipticite dynamique
      Write(*,*) 'Relative change of dynamical ellipticity'
      Write(*,*) 'ENTER gamma/gamma_0  (default 1) '
      read (*,*) FGAM
      write (*,*) 'FGAM = ', FGAM
      CP1 = CP1*FGAM
      CP2 = CP2*FGAM
      CP3 = CP3*FGAM**2
      CP4 = CP4*FGAM
      write (*,*)
***-------------------- effet de maree
      Write(*,*) 'Relative tidal effet'
      write (*,*) ' 0 : no tidal effect '
      write (*,*) ' 1 :  -4.6 D-18 seconds**-1'
      Write(*,*) 'ENTER cmar  (default 0) '
      read (*,*) CMAR
      write (*,*) 'CMAR = ', CMAR
      write (*,*)
C CONSTANTES POUR PRENDRE EN COMPTE L'EFFET DE MAREES
      AK2=-4.6D-18*86400D0*365.25D0*CMAR
      AK1=51D0*AK2*AML/OM
***-----------------------------------------------------------
      IF (IPT.EQ.1) THEN
         write (*,*) 
         write (*,*) 'Informations for internal check'
         write (*,*)
         PRINT *,'VITESSE ANGULAIRE DE LA TERRE',OM*RDSC
         PRINT *,'PRECESSION EN ASCENSION DROITE',APP*RDSC*100
         PRINT *,'ELD',ELD
         PRINT *,'RFL0',RFL0
         PRINT *,'RFL1',RFL1
         PRINT *,'RFL3',RFL3
         PRINT *,'RFS',RFS
         PRINT *,'RFL0+RFL1+RFL3+RFS',RFL0+RFL1+RFL3+RFS
         PRINT *,'CP1,CP2,CP3,CP4  EN SECONDES PAR AN '
         PRINT *,CP1*RDSC,CP2*RDSC,CP3*RDSC,CP4*RDSC
         PRINT *,'AK1',AK1
         PRINT *,'AK2',AK2
         write (*,*)
      ENDIF
      RETURN
      END
      SUBROUTINE PRECES(t,AK,AH,AQ,AP,DK,DH,DQ,DP,AKI,DKI)
***********************************************************************
*   SECOND MEMBRE DES EQUATIONS DE LA PRECESSION POUR LA TERRE        *
*   CORRIGE POUR TENIR COMPTE DES EFFETS DE MAREES                    *
*                                                                     *
*  (c) Astronomie et Systemes Dynamiques, Bureau des Longitudes (1993)*
***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Y)
      DIMENSION AKI(2),DKI(2)
      COMMON/CONST/APG,CP1,CP2,CP3,CP4,AK1,AK2
      AKIX=AKI(1)
      AKIY=AKI(2)
C  X=COS(EPS)
      X=DSQRT(1.D0-(AKIX**2 + AKIY**2))
      SINEPS=DSQRT(AKIX**2 + AKIY**2)
      RS=(1-AK**2-AH**2)**(-1.5D0)/2-0.522D-6
      COR1=1.D0+AK1*T
      COR2=1.D0+2.D0*AK2*T
C
      RTOT=COR1*(COR2*(CP1*X+CP2*(2.D0*X**2-1.D0)/SINEPS
     $     +COR2*CP3*(6.D0*X**2-1.D0))+CP4*RS*X)
C
      CC=AQ*DP-AP*DQ
      DD=2/DSQRT(1-AP**2-AQ**2)
      AA=DD*(DQ+AP*CC)
      BB=DD*(DP-AQ*CC)
      COEF=RTOT-2.D0*CC-APG
      DKIY=AKIX*COEF-BB*X
      DKIX=-AKIY*COEF+AA*X
      DKI(1)=DKIX
      DKI(2)=DKIY
      END