subroutine calcmilank(tim,obliq,perid,eccen)
          
      use planet_h          

      implicit none

#include "comcstfi.h"

!==================================================================
!     Purpose
!     -------
!     calculate the milankovitch parametersdepeding on the
!     date (tim, in Myr)
!     Obliquity
!     Day of perihelion from the vernal equinox 
!     Eccentricity

!     Inputs
!     ------ 
!     tim                   Date (Myr)
!     
!     Outputs
!     -------
!     obliq                 new value for obliquity 
!     perid                 new value for peri_day 
!     eccen                 new value for eccentricity
!     
!     Authors
!     ------- 
!     Tanguy Bertrand (new equations in 2022)
!     
!==================================================================

!------------------------------------------------------------------
!     Arguments

      REAL,intent(in) :: tim
      REAL,intent(out) :: perid,obliq,eccen
      real teta,xo,xref,psi,phi
!------------------------------------------------------------------

      write(*,*) 'TB22 tim  =',tim

      psi=72.8+91.*tim
      phi=19.5+130.2*tim

      ! update new obliquity
      obliq=115.5+11.8*sin(phi*pi/180.)

      ! update new eccentricity
      eccen=0.244+0.022*cos(psi*pi/180.)+0.005*cos(3*psi*pi/180.)

      ! update new solar longitude of perihelion
      teta=pi/180.*modulo(19.5+130.2*tim-24.*sin(psi*pi/180.),360.)
      xo=2.*atan(((1-eccen)/(1+eccen))**0.5*tan(teta/2.))
      xref=xo-eccen*sin(xo)
      perid=year_day*(1.-xref/(2.*pi))

      if (perid.ge.year_day) then
         perid = perid - year_day
      endif
      write(*,*) 'new ecc, lsp, perid=',eccen,360.-teta*180./pi,perid
      write(*,*) 'new obliquity=',obliq

      end subroutine calcmilank