MODULE recomp_orb_param_mod
!=======================================================================
!
! Purpose: Recompute orbit parameters based on values by Laskar et al.,
!
! Author: RV 
!=======================================================================
      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 :: timeperi_ls	        ! time of peri in ls
      integer ilask			! Loop variable
      real :: halfaxe                   ! Million of km
      real :: unitastr

      real xa,xb,xc,ya,yb


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

      aphelie=apoastr
      periheli=periastr
#endif

          Year=year_bp_ini+year_PEM+final_iter

          timeperi_ls=timeperi*360/2/pi

!          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=", timeperi_ls
          print *, "recomp_orb_param, Old timeperi=", timeperi

        xc=Year

        do ilask=last_ilask,1,-1

          xa=yearlask(ilask)
          xb=yearlask(ilask+1)

           if(yearlask(ilask) .GT.Year) then
             if(var_obl) then
              ya=oblask(ilask)
              yb=oblask(ilask+1)
              obliquit=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
             endif
             if(var_ex) then
              ya=exlask(ilask)
              yb=exlask(ilask+1)
              e_elips=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
             endif
             if(var_lsp) then
               if(lsplask(ilask)-lsplask(ilask+1) .gt.200) then
                 if(lsplask(ilask).lt.lsplask(ilask+1)) then
                   lsplask(ilask)=lsplask(ilask)+360
                 else
                   lsplask(ilask+1)=lsplask(ilask+1)+360
                 endif
               endif

              ya=lsplask(ilask)
              yb=lsplask(ilask+1)
              timeperi_ls=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
              timeperi_ls=modulo(timeperi_ls,360.)
              halfaxe=227.94
              timeperi=timeperi_ls*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
             endif
              exit
           endif
        enddo

       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