MODULE orbit_param_criterion_mod

!=======================================================================
!
! Purpose: Compute the maximum number of iteration for PEM based on the 
! stopping criterion given by the orbital parameters
!
! Author: RV, JBC
!=======================================================================

IMPLICIT NONE

CONTAINS

SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)

#ifdef CPP_IOIPSL
    use IOIPSL,          only: getin
#else
    ! if not using IOIPSL, we still need to use (a local version of) getin
    use ioipsl_getincom, only: getin
#endif
#ifndef CPP_STD
    use planete_h,   only: e_elips, obliquit, lsperi
#else
    use planete_mod, only: e_elips, obliquit, lsperi
#endif
  use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
  use comconst_mod,   only: pi
  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask

  IMPLICIT NONE

!--------------------------------------------------------
! Input Variables
!--------------------------------------------------------
  integer, intent(in) :: i_myear ! Martian year of the simulation

!--------------------------------------------------------
! Output Variables
!--------------------------------------------------------
  integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much

!--------------------------------------------------------
! Local variables 
!--------------------------------------------------------
  integer :: Year                                           ! Year of the simulation
  real    :: Lsp                                            ! Ls perihelion
  real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
  real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
  real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
  integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
  integer :: ilask, iylask                                  ! Loop variables
  real    :: xa, xb, ya, yb, yc_obl, yc_ecc, yc_lsp         ! Variables for interpolation

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

call ini_lask_param_mod ! Allocation of variables

write(*,*) 'orbit_param_criterion, Current year =', Year
write(*,*) 'orbit_param_criterion, Current obl. =', obliquit
write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips
write(*,*) 'orbit_param_criterion, Current Lsp  =', Lsp

! We read the file
open(73,file = 'obl_ecc_lsp.asc')
do ilask = 1,size(yearlask,1)
    read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask)
    yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
    if (yearlask(ilask) > Year) last_ilask = ilask + 1
enddo
close(73)
write(*,*) 'Corresponding line in the obl_ecc_lsp.asc file =', last_ilask

! Constant max change case 
max_change_obl = 0.6
max_change_ecc = 2.8e-3
max_change_lsp = 18.

call getin('max_change_obl', max_change_obl)
max_obl = obliquit + max_change_obl
min_obl = obliquit - max_change_obl
write(*,*) 'Maximum obliquity accepted =', max_obl
write(*,*) 'Minimum obliquity accepted =', min_obl

call getin('max_change_ecc', max_change_ecc)
max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ex < 1.
min_ecc = max(e_elips - max_change_ecc,0.) ! ex >= 0.
write(*,*) 'Maximum eccentricity accepted =', max_ecc
write(*,*) 'Minimum eccentricity accepted =', min_ecc

call getin('max_change_lsp', max_change_lsp)
max_lsp = Lsp + max_change_lsp
min_lsp = Lsp - max_change_lsp
write(*,*) 'Maximum Lsp accepted =', max_lsp
write(*,*) 'Minimum Lsp accepted =', min_lsp
! End Constant max change case 

! If we do not want some orb parameter to change, they should not be a stopping criterion,
! So the number of iteration corresponding is set to maximum
max_obl_iter = 1000000
max_ecc_iter = 1000000
max_lsp_iter = 1000000

        ! Tendency of the orbital parameter for the considered year gives the limitation between min and max
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)
            if (ya < yb) then ! Increasing -> max is the limitation
                yc_obl = max_obl
            else ! Decreasing -> min is the limitation
                yc_obl = min_obl
            endif
        endif

        ! Eccentricity
        if (var_ecc) then
            ya = ecclask(ilask)
            yb = ecclask(ilask - 1)
            if (ya < yb) then ! Increasing -> max is the limitation
                yc_ecc = max_ecc
            else ! Decreasing -> min is the limitation
                yc_ecc = min_ecc
            endif
        endif

        ! Lsp
        if (var_lsp) then
            ya = lsplask(ilask)
            yb = lsplask(ilask - 1)
             if (ya < yb) then ! Increasing -> max is the limitation
                if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around
                    yc_lsp = min_lsp
                else
                    yc_lsp = max_lsp
                endif
            else ! Decreasing -> min is the limitation
                if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around
                    yc_lsp = max_lsp
                else
                    yc_lsp = min_lsp
                endif
            endif
        endif
        iylask = ilask
        exit ! The loop is left as soon as the right interval is found
    endif
enddo

if (ilask == 1) then
    write(*,*) 'The year does not match with Laskar data in the file obl_ecc_lsp.asc.'
    stop
endif

! Linear interpolation gives the maximum reachable year according to the limitation
! Obliquity
if (var_obl) then
    do ilask = iylask,2,-1
        ya = obllask(ilask)
        yb = obllask(ilask - 1)
        if (min(ya,yb) <= yc_obl .and. yc_obl < max(ya,yb)) then
            xa = yearlask(ilask)
            xb = yearlask(ilask - 1)
            max_obl_iter = floor((yc_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
            exit ! The loop is left as soon as the right interval is found
        endif
    enddo
endif
! Eccentricity
if (var_ecc) then
    do ilask = iylask,2,-1
        ya = ecclask(ilask)
        yb = ecclask(ilask - 1)
        if (min(ya,yb) <= yc_ecc .and. yc_ecc < max(ya,yb)) then
            xa = yearlask(ilask)
            xb = yearlask(ilask - 1)
            max_ecc_iter = floor((yc_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
            exit ! The loop is left as soon as the right interval is found
        endif
    enddo
endif
! Lsp
if (var_lsp) then
    do ilask = iylask,2,-1
        ya = lsplask(ilask)
        yb = lsplask(ilask - 1)
        if (min(ya,yb) <= yc_lsp .and. yc_lsp < max(ya,yb)) then
            xa = yearlask(ilask)
            xb = yearlask(ilask - 1)
            max_lsp_iter = floor((yc_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
            exit ! The loop is left as soon as the right interval is found
        endif
    enddo
endif

write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter
write(*,*) 'Max. number of iterations for the ecc. parameter =', max_ecc_iter
write(*,*) 'Max. number of iterations for the Lsp parameter  =', max_lsp_iter

year_iter_max = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
if (year_iter_max > 0) then
    write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
else
    write(*,*) 'The max. number of iterations is not compatible with Laskar data in the file obl_ecc_lsp.asc.'
    stop
endif

END SUBROUTINE orbit_param_criterion

!********************************************************************************   

END MODULE orbit_param_criterion_mod