MODULE recomp_orb_param_mod
!=======================================================================
!
! Purpose: Recompute orbit parameters based on values by Laskar et al.,
!
! Author: RV, JBC
!=======================================================================
      IMPLICIT NONE

      CONTAINS

      SUBROUTINE recomp_orb_param(final_iter)

      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
#ifndef CPP_STD      
      USE comconst_mod, ONLY: pi
      USE planete_h, ONLY: e_elips, obliquit, timeperi, periheli,aphelie,p_elips,peri_day,year_day
#else
      use planete_mod, only: e_elips, obliquit, timeperi,periastr,apoastr,p_elips,peri_day,year_day
      USE comcstfi_mod, only: pi

#endif

      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
   				end_lask_param_mod,last_ilask

      IMPLICIT NONE

!--------------------------------------------------------
! Input Variables
!--------------------------------------------------------

      integer,intent(in) :: final_iter ! Number of year iteration of the PEM 

!--------------------------------------------------------
! Output Variables
!--------------------------------------------------------

!--------------------------------------------------------
! Local variables 
!--------------------------------------------------------

      real :: Year     ! Year of the simulation
      real :: Lsp      ! time of peri in ls
      integer ilask    ! Loop variable
      real :: halfaxe  ! Million of km
      real :: unitastr

      real xa,xb,ya,yb


      ! **********************************************************************
      ! 0. Initializations
      ! **********************************************************************
#ifdef CPP_STD 
      real aphelie
      real periheli

      aphelie=apoastr
      periheli=periastr
#endif

          Year = year_bp_ini + year_PEM + final_iter ! We compute the current year
          Lsp = 360. - timeperi*360./(2.*pi)         ! We convert in degree

!          Year=-953200
          print*, "recomp_orb_param, Old year=", year_bp_ini+year_PEM
          print*, "recomp_orb_param, New year=", Year
          print*, "recomp_orb_param, Old obl=", obliquit
          print*, "recomp_orb_param, Old ex=", e_elips
          print*, "recomp_orb_param, Old lsp=", Lsp
          print*, "recomp_orb_param, Old timeperi=", timeperi

        do ilask = last_ilask,2,-1
          xa = yearlask(ilask)
          xb = yearlask(ilask - 1)
          if (xa <= Year .and. Year < xb) then
            ! Obliquity
            if (var_obl) then
              ya = oblask(ilask)
              yb = oblask(ilask - 1)
              obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
            endif

            ! Excentricity
            if (var_ex) then
              ya = exlask(ilask)
              yb = exlask(ilask - 1)
              e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
            endif

            ! Lsp
            if (var_lsp) then
              ya = lsplask(ilask)
              yb = lsplask(ilask - 1)
              if (yb - ya > 300.) then ! If modulo is "crossed"
                if (ya < yb) then ! Increasing
                  yb = yb - 360.
                else ! Decreasing
                  yb = yb + 360.
                endif
              endif
              Lsp = (Year - xa)*(yb - ya)/(xb - xa) + ya
              Lsp = modulo(Lsp,360.)
            endif
            exit ! The loop is left as soon as the right interval is found
          endif
        enddo
        halfaxe = 227.94
        timeperi = Lsp*2.*pi/360.
        periheli = halfaxe*(1 - e_elips)
        aphelie =  halfaxe*(1 + e_elips)
        unitastr = 149.597927
        p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
#ifndef CPP_STD  
        call call_dayperi(timeperi,e_elips,peri_day,year_day)
#endif

       print*, "recomp_orb_param, Final year of the PEM run:", Year
       print*, "recomp_orb_param, New obl=", obliquit
       print*, "recomp_orb_param, New ex=", e_elips
       print*, "recomp_orb_param, New timeperi=", timeperi

      END SUBROUTINE recomp_orb_param

!********************************************************************************   
      
      END MODULE recomp_orb_param_mod