source: trunk/LMDZ.COMMON/libf/evolution/orbit.F90 @ 4003

MODULE orbit
!-----------------------------------------------------------------------
! NAME
!     orbit
! DESCRIPTION
!     Compute orbital-parameter-based stopping criteria and update
!     planetary orbit parameters.
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! MODULE VARIABLES
! ----------------
real, dimension(:), allocatable :: year_lask  ! year before present from Laskar et al.
real, dimension(:), allocatable :: obl_lask   ! obliquity [deg]
real, dimension(:), allocatable :: ecc_lask   ! eccentricity
real, dimension(:), allocatable :: lsp_lask   ! ls perihelie [deg]
integer                         :: iyear_lask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE read_orbitpm(Year)
!-----------------------------------------------------------------------
! NAME
!     read_orbitpm
!
! DESCRIPTION
!     Read Laskar orbital file and prepare arrays; locate index for
!     closest lower year relative to current simulation year.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!     Converts kilo Earth years to Martian years using 'convert_years'.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use evolution, only: convert_years

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real, intent(in) :: Year ! Martian year of the simulation

! LOCAL VARIABLES
! ---------------
integer :: n ! Number of lines in Laskar files
integer :: ierr, i

! CODE
! ----
n = 0
open(1,file = 'obl_ecc_lsp.asc',status = 'old',action = 'read')
do
    read(1,*,iostat = ierr)
    if (ierr /= 0) exit
    n = n + 1
enddo
close(1)
allocate(year_lask(n),obl_lask(n),ecc_lask(n),lsp_lask(n))

iyear_lask = 0
do i = 1,size(year_lask,1)
    read(1,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i)
    year_lask(i) = year_lask(i)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
    if (year_lask(i) > Year) iyear_lask = i + 1
enddo
close(1)
if (iyear_lask == 0 .or. iyear_lask == size(year_lask,1) + 1) then
    error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".'
else
    write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', iyear_lask
endif

END SUBROUTINE read_orbitpm
!=======================================================================

!=======================================================================
SUBROUTINE end_orbit
!-----------------------------------------------------------------------
! NAME
!     end_orbit
!
! DESCRIPTION
!     Deallocate orbital data arrays.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! CODE
! ----
if (allocated(year_lask)) deallocate(year_lask)
if (allocated(obl_lask))  deallocate(obl_lask)
if (allocated(ecc_lask))  deallocate(ecc_lask)
if (allocated(lsp_lask))  deallocate(lsp_lask)

END SUBROUTINE end_orbit
!=======================================================================

!=======================================================================
SUBROUTINE orbit_param_criterion(i_myear,nyears_maxorb)
!-----------------------------------------------------------------------
! NAME
!     orbit_param_criterion
!
! DESCRIPTION
!     Determine maximum PEM iterations before orbital params change
!     beyond admissible thresholds.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------
#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 planete_h,      only: e_elips, obliquit, lsperi
use phys_constants, only: pi
use evolution,      only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real, intent(in)  :: i_myear       ! Martian year of the simulation
real, intent(out) :: nyears_maxorb ! Maximum number of iteration of the PEM before orb changes too much

! LOCAL VARIABLES
! ---------------
real                            :: 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
real                            :: nyears_maxobl, nyears_maxecc, nyears_maxlsp    ! Maximum year iteration before reaching an inadmissible value of orbit param
real                            :: xa, xb, ya, yb                  ! Variables for interpolation
logical                         :: found_obl, found_ecc, found_lsp ! Flag variables for orbital parameters
real                            :: lsp_corr                        ! Lsp corrected if the "modulo is crossed"
real                            :: delta                           ! Intermediate variable
real, dimension(:), allocatable :: lsplask_unwrap                  ! Unwrapped sequence of Lsp from Laskar's file
integer                         :: i

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

call read_orbitpm(Year) ! 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

! Unwrap the Lsp changes in Laskar's file
allocate(lsplask_unwrap(iyear_lask))
lsplask_unwrap(iyear_lask) = lsp_lask(iyear_lask)
do i = iyear_lask,2,-1
    delta = lsp_lask(i - 1) - lsp_lask(i)
    if (delta > 180.) then
        lsp_corr = lsp_lask(i - 1) - 360.
    else if (delta < -180.) then
        lsp_corr = lsp_lask(i - 1) + 360.
    else
        lsp_corr = lsp_lask(i - 1)
    endif
    lsplask_unwrap(i - 1) = lsplask_unwrap(i) + lsp_corr - lsp_lask(i)
enddo

! Correct the current Lsp according to the unwrapping
delta = Lsp - lsplask_unwrap(iyear_lask)
if (delta > 180.) then
    Lsp = Lsp - 360.
else if (delta < -180.) then
    Lsp = Lsp + 360.
endif

! Constant max change case
max_change_obl = 1.
max_change_ecc = 5.e-3
max_change_lsp = 20.

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) ! ecc < 1.
min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 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
nyears_maxobl = 1000000.
nyears_maxecc = 1000000.
nyears_maxlsp = 1000000.

! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
found_obl = .false.
found_ecc = .false.
found_lsp = .false.
if (.not. var_obl) found_obl = .true.
if (.not. var_ecc) found_ecc = .true.
if (.not. var_lsp) found_lsp = .true.

do i = iyear_lask,2,-1
    xa = year_lask(i)
    xb = year_lask(i - 1)

    ! Obliquity
    if (var_obl .and. .not. found_obl) then
        ya = obl_lask(i)
        yb = obl_lask(i - 1)
        if (yb < min_obl) then ! The minimum accepted is overtaken
            nyears_maxobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_obl = .true.
        else if (max_obl < yb) then ! The maximum accepted is overtaken
            nyears_maxobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_obl = .true.
        endif
    endif

    ! Eccentricity
    if (var_ecc .and. .not. found_ecc) then
        ya = ecc_lask(i)
        yb = ecc_lask(i - 1)
        if (yb < min_ecc) then ! The minimum accepted is overtaken
            nyears_maxecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_ecc = .true.
        else if (max_ecc < yb) then ! The maximum accepted is overtaken
            nyears_maxecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_ecc = .true.
        endif
    endif

    ! Lsp
    if (var_lsp .and. .not. found_lsp) then
        ya = lsplask_unwrap(i)
        yb = lsplask_unwrap(i - 1)
        if (yb < min_lsp) then ! The minimum accepted is overtaken
            nyears_maxlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_lsp = .true.
        else if (max_lsp < yb) then ! The maximum accepted is overtaken
            nyears_maxlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_lsp = .true.
        endif
    endif

    if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found
enddo

deallocate(lsplask_unwrap)

if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".'
if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".'
if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".'

write(*,*) 'Max. number of years for the obl. parameter =', nyears_maxobl
write(*,*) 'Max. number of years for the ecc. parameter =', nyears_maxecc
write(*,*) 'Max. number of years for the Lsp parameter  =', nyears_maxlsp

nyears_maxorb = min(nyears_maxobl,nyears_maxecc,nyears_maxlsp)
if (nyears_maxorb < 1.) nyears_maxorb = 1. ! nyears_maxorb should be at least equal to 1
write(*,*) 'So the max. number of years for the orbital criteria =', nyears_maxorb

END SUBROUTINE orbit_param_criterion
!=======================================================================

!=======================================================================
SUBROUTINE set_new_orbitpm(i_myear,i_myear_leg)
!-----------------------------------------------------------------------
! NAME
!     set_new_orbitpm
!
! DESCRIPTION
!     Recompute orbit parameters from Laskar values at a new year and
!     update planetary constants; handles Lsp modulo crossings.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

use evolution,        only: year_bp_ini, var_obl, var_ecc, var_lsp
use phys_constants,   only: pi
use planete_h,        only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
use call_dayperi_mod, only: call_dayperi

! DECLARATION
! -----------
implicit none

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

! LOCAL VARIABLES
! ---------------
real    :: Year           ! Year of the simulation
real    :: Lsp            ! Ls perihelion
real    :: halfaxe        ! Million of km
integer :: i
real    :: unitastr
real    :: xa, xb, ya, yb ! Variables for interpolation
logical :: found_year     ! Flag variable

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

write(*,*) 'set_new_orbitpm, Old year =', Year - i_myear_leg
write(*,*) 'set_new_orbitpm, Old obl. =', obliquit
write(*,*) 'set_new_orbitpm, Old ecc. =', e_elips
write(*,*) 'set_new_orbitpm, Old Lsp  =', Lsp
write(*,*) 'set_new_orbitpm, New year =', Year

found_year = .false.
if (Year < 0.) then
    do i = iyear_lask,2,-1
        xa = year_lask(i)
        xb = year_lask(i - 1)
        if (xa <= Year .and. Year < xb) then
            ! Obliquity
            if (var_obl) then
                ya = obl_lask(i)
                yb = obl_lask(i - 1)
                obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
            endif

            ! Eccentricity
            if (var_ecc) then
                ya = ecc_lask(i)
                yb = ecc_lask(i - 1)
                e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
            endif

            ! Lsp
            if (var_lsp) then
                ya = lsp_lask(i)
                yb = lsp_lask(i - 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
else if (Year >= 0 .and. Year < 100.) then ! For present orbital characteristics
    if (var_obl) obliquit = 25.19
    if (var_ecc) e_elips = 0.0934
    if (var_lsp) Lsp = 251.
    found_year = .true.
endif

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

write(*,*) 'set_new_orbitpm, New obl. =', obliquit
write(*,*) 'set_new_orbitpm, New ecc. =', e_elips
write(*,*) 'set_new_orbitpm, New Lsp  =', Lsp

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

call call_dayperi(Lsp,e_elips,peri_day,year_day)

call end_orbit()

END SUBROUTINE set_new_orbitpm
!=======================================================================

END MODULE orbit