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(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

      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
#ifndef CPP_STD
      USE planete_h, ONLY: e_elips, obliquit, timeperi
#else
      use planete_mod, only: e_elips, obliquit, timeperi
#endif
      USE comconst_mod, ONLY: pi
      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
   				ini_lask_param_mod,last_ilask

      IMPLICIT NONE

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

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

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

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

      real :: Year	                ! Year of the simulation
      real :: Lsp	                ! Year of the simulation
      integer nlask, ilask, iylask      ! Loop variables
      parameter (nlask = 20001)         ! Number of line in Laskar file

      real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible

      real max_obl,max_ex,max_lsp       ! Maximum value of orbit param given the acceptable percentage 
      real min_obl,min_ex,min_lsp       ! Maximum value of orbit param given the acceptable percentage
      real max_obl_iter,max_ex_iter,max_lsp_iter ! Maximum year iteration before reaching an unacceptable value
      real xa,xb,ya,yb,yc

      ! **********************************************************************
      ! 0. Initializations
      ! **********************************************************************

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

	  call ini_lask_param_mod(nlask) ! Allocation of variables

          print*, "orbit_param_criterion, Year in pem.def=", year_bp_ini
          print*, "orbit_param_criterion, Year in the startpem.nc =", year_PEM
          print*, "orbit_param_criterion, Current year=", Year
          print*, "orbit_param_criterion, Current obl=", obliquit
          print*, "orbit_param_criterion, Current ex=", e_elips
          print*, "orbit_param_criterion, Current lsp=", Lsp

! We read the file

          open(73,file='ob_ex_lsp.asc')
          do ilask = 1,nlask
            read(73,*) yearlask(ilask), oblask(ilask), exlask(ilask), lsplask(ilask)
            yearlask(ilask) = yearlask(ilask)*1000.
            if (yearlask(ilask) > Year) last_ilask = ilask + 1
          enddo
          close(73)

       print*, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask

!Constant max change case 

        max_change_obl = 0.5
        max_change_ex = 0.1
        max_change_lsp = 40.

        CALL getin('max_change_obl', max_change_obl)
	max_obl = obliquit + max_change_obl
	min_obl = obliquit - max_change_obl
        print*, "Maximum obliquity accepted=", max_obl
        print*, "Minimum obliquity accepted=", min_obl

        CALL getin('max_change_ex', max_change_ex)
	max_ex = e_elips + max_change_ex
	min_ex = e_elips - max_change_ex
        print*, "Maximum excentricity accepted=", max_ex
        print*, "Minimum excentricity accepted=", min_ex

        CALL getin('max_change_lsp', max_change_lsp)
	max_lsp = Lsp + max_change_lsp
	min_lsp = Lsp - max_change_lsp
        print*, "Maximum lsp accepted=", max_lsp
        print*, "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 = 999999
        max_ex_iter = 999999
        max_lsp_iter = 999999

        ! 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 = oblask(ilask)
              yb = oblask(ilask - 1)
              if (ya < yb) then ! Increasing -> max is the limitation
                yc = max_obl
              else ! Decreasing -> min is the limitation
                yc = min_obl
              endif
            endif

            ! Excentricity
            if (var_ex) then
              ya = exlask(ilask)
              yb = exlask(ilask - 1)
              if (ya < yb) then ! Increasing -> max is the limitation
                yc = max_ex
              else ! Decreasing -> min is the limitation
                yc = min_ex
              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
                  !yb = yb - 360.
                  yc = min_lsp
                else
                  yc = 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
                  !yb = yb + 360.
                  yc = max_lsp
                else
                  yc = 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 ob_ex_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 = oblask(ilask)
            yb = oblask(ilask - 1)
            if (ya <= yc .and. yc < yb) then
              xa = yearlask(ilask)
              xb = yearlask(ilask - 1)
              max_obl_iter = int((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
              exit ! The loop is left as soon as the right interval is found
            endif
          enddo
	endif
        ! Excentricity
        if (var_ex) then
          do ilask = iylask,2,-1
            ya = exlask(ilask)
            yb = exlask(ilask - 1)
            if (ya <= yc .and. yc < yb) then
              xa = yearlask(ilask)
              xb = yearlask(ilask - 1)
              max_ex_iter = int((max_ex - 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 (ya <= yc .and. yc < yb) then
              xa = yearlask(ilask)
              xb = yearlask(ilask - 1)
              max_lsp_iter = int((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
              exit ! The loop is left as soon as the right interval is found
            endif
          enddo
	endif

      print*, "Maximum number of iteration for the obl. parameter=", max_obl_iter
      print*, "Maximum number of iteration for the ex. parameter=", max_ex_iter
      print*, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter

      year_iter_max = min(max_obl_iter,max_ex_iter,max_lsp_iter)
      if (year_iter_max > 0) then
        print*, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max
      else
        write(*,*) 'The max. number of iteration (year) is not compatible with Laskar data in the file ob_ex_lsp.asc.'
        stop
      endif

        END SUBROUTINE orbit_param_criterion

!********************************************************************************   
      
      END MODULE orbit_param_criterion_mod