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(i_myear,year_iter)

  use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
#ifndef CPP_STD      
    use comconst_mod, only: pi
    use planete_h,    only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
#else
    use planete_mod,  only: e_elips, obliquit, lsperi, periastr, apoastr, p_elips, peri_day, year_day
    use comcstfi_mod, only: pi
#endif

  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask

  IMPLICIT NONE

!--------------------------------------------------------
! Input Variables
!--------------------------------------------------------
  integer, intent(in) :: i_myear   ! Number of simulated Martian years
  integer, intent(in) :: year_iter ! Number of iterations of the PEM

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

!--------------------------------------------------------
! Local variables 
!--------------------------------------------------------
  integer :: Year           ! Year of the simulation
  real    :: Lsp            ! Ls perihelion
  integer :: ilask          ! Loop variable
  real    :: halfaxe        ! Million of km
  real    :: unitastr
  real    :: xa, xb, ya, yb ! Variables for interpolation

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

    aphelie=apoastr
    periheli=periastr
#endif

Year = year_bp_ini + i_myear ! We compute the current year
Lsp = lsperi*180./pi         ! We convert in degree

write(*,*) 'recomp_orb_param, Old year =', Year - year_iter
write(*,*) 'recomp_orb_param, Old obl. =', obliquit
write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp

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 = obllask(ilask)
            yb = obllask(ilask - 1)
            obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
        endif

        ! Eccentricity
        if (var_ecc) then
            ya = ecclask(ilask)
            yb = ecclask(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
Lsp = Lsp*pi/180.
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(Lsp,e_elips,peri_day,year_day)
#endif

write(*,*) 'recomp_orb_param, New year =', Year
write(*,*) 'recomp_orb_param, New obl. =', obliquit
write(*,*) 'recomp_orb_param, New ecc. =', e_elips
write(*,*) 'recomp_orb_param, New Lsp  =', Lsp

END SUBROUTINE recomp_orb_param

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