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,i_myear_leg)

use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp
use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
#ifndef CPP_STD
    use comcstfi_h,   only: pi
    use planete_h,    only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
#else
    use comcstfi_mod, only: pi
    use planete_mod,  only: e_elips, obliquit, lsperi, periastr, apoastr, p_elips, peri_day, year_day
#endif

implicit none

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

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

!--------------------------------------------------------
! Local variables
!--------------------------------------------------------
real    :: 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
logical :: found_year     ! Flag variable

! **********************************************************************
! 0. Initializations
! **********************************************************************
#ifdef CPP_STD
    real :: aphelie, 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 - i_myear_leg
write(*,*) 'recomp_orb_param, Old obl. =', obliquit
write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp

found_year = .false.
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 (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval
                if (ya < yb) then ! Lsp should be decreasing
                    ya = ya + 360.
                else ! Lsp should be increasing
                    yb = yb + 360.
                endif
            endif
            Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.)
        endif
        found_year = .true.
        exit ! The loop is left as soon as the right interval is found
    endif
enddo

if (.not. found_year) error stop 'The new year could not be found in the file "obl_ecc_lsp.asc".'

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