program climavar
*--------------------------------------------------------------------------
*     This small program convert the ASCII files 
*     of the orbital elements          t, k, h, q, p
*     and of the precession quantities t, eps, phi
*
*     into ASCII files for the usual climatic variables 
* 
*     t, e, eps ,e*sin(pibar)
*
*     where eps is the obliquity, e the eccentricity, 
*     and   pibar the longitude of perihelion from moving equinox
*     e*sin(pibar) is thus the climatic precession
* 
*     this program can be use as an example to derive other 
*     quantities from the existing ASCII
*
*
*      parameters are taken in the  NAMELIST file : 'climavar.par'
*
*       nomascpos   :  ASCII file for elements t,k,h,q,p
*                      positive time
*       nomascneg   :  ASCII file for elements t,h,h,q,p
*                      negative time
*       nomprecpos  :  ASCII file for elements t,eps,phi
*                      positive time
*       nomprecneg  :  ASCII file for elements t,eps,phi
*                      negative time
*       datedebut   :  starting time in millions of years
*       datefin     :  end      time in millions of years
*                        ( -20 <= datedebut < datefin <= +10 )
*       nomsolpos   :  ASCII file for elements t,e,eps,e*sin(pibar)
*                      positive time
*       nomsolneg   :  ASCII file for elements t,e,eps,e*sin(pibar)
*                      negative time
*
*     (c) ASD/BdL version 0.8 (4/6/93)
*-------------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      parameter(nbel=5)
      character *20 statut
      character *50 nomascpos,nomascneg,nomprecpos,nomprecneg,
     &nomsolpos,nomsolneg
      NAMELIST/NAMSTD/nomascpos,nomascneg,nomprecpos,nomprecneg,
     &nomsolpos,nomsolneg,datedebut,datefin,statut
      open(8,file='climavar.par',form='formatted',status='old')
      read(8,NAMSTD)
      close(8)
      pas=1D3
      if (datedebut.lt.0.D0) then
        open(11,file=nomascneg,status='old')
        open(12,file=nomprecneg,status='old')
        open(13,file=nomsolneg,status=statut)
        nblignes=int(abs(datedebut*1D6/pas))+1
        do i=1,nblignes
           read(11,*) t,rk,rh,rq,rp
           read(12,*) t,eps,rphi
           e=sqrt(rh**2+rk**2)
           pitilde=atan2(rh,rk)
           pibar=pitilde+rphi
           write(13,1000) t,e,eps,e*sin(pibar)
        enddo
        close(11)
        close(12)
        close(13)
      endif
      if (datefin.gt.0.D0) then
        open(11,file=nomascpos,status='old')
        open(12,file=nomprecpos,status='old')
        open(13,file=nomsolpos,status=statut)
        nblignes=int(abs(datefin*1D6/pas))+1
        do i=1,nblignes
           read(11,*) t,rk,rh,rq,rp
           read(12,*) t,eps,rphi
           e=sqrt(rh**2+rk**2)
           pitilde=atan2(rh,rk)
           pibar=pitilde+rphi
           write(13,1000) t,e,eps,e*sin(pibar)
        enddo
        close(11)
        close(12)
        close(13)
      endif
1000  format(1x,f10.0,3d20.12)
      end