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

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

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4, largest_nb, minieps

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

! PARAMETERS
! ----------
! File-related:
character(15),                       parameter :: orbitfile_name = 'obl_ecc_lsp.asc'
integer(di),                         protected :: iyear_lask! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
real(dp), dimension(:), allocatable, protected :: year_lask ! Current_year before present from Laskar et al.
real(dp), dimension(:), allocatable, protected :: obl_lask  ! Obliquity in Laskar's data [deg]
real(dp), dimension(:), allocatable, protected :: ecc_lask  ! Eccentricity in Laskar's data
real(dp), dimension(:), allocatable, protected :: lsp_lask  ! Ls perihelie in Laskar's data [deg]
! Planet-related:
real(dp),                            protected :: obliquity    ! Planet obliquity [deg]
real(dp),                            protected :: eccentricity ! Orbit eccentricity
real(dp),                            protected :: ls_peri      ! Solar longitude of the perihelion [deg]
real(dp),                            protected :: perihelion   ! Perihelion [Mkm]
real(dp),                            protected :: semimaj_axis ! Semi-major axis [Mkm]
real(dp),                            protected :: aphelion     ! Aphelion [Mkm]
real(dp),                            protected :: date_peri    ! Date of perihelion [sol]
real(dp),                            protected :: sol_len_s    ! Length of a sol [s]
real(dp),                            protected :: yr_len_sol   ! Length of a year [sol]
! Evolution-related:
logical(k4), protected :: evol_orbit     ! Flag to follow the orbital parameters
logical(k4), protected :: evol_obl       ! Flag the PEM to follow obliquity
logical(k4), protected :: evol_ecc       ! Flag the PEM to follow eccenticity
logical(k4), protected :: evol_lsp       ! Flag the PEM to follow Ls perihelie
real(dp),    protected :: max_change_obl ! Maximum admissible change for obliquity
real(dp),    protected :: max_change_ecc ! Maximum admissible change for eccentricity
real(dp),    protected :: max_change_lsp ! Maximum admissible change for Lsp

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

!=======================================================================
SUBROUTINE set_orbit_config(evol_orbit_in,evol_obl_in,evol_ecc_in,evol_lsp_in,max_change_obl_in,max_change_ecc_in,max_change_lsp_in)
!-----------------------------------------------------------------------
! NAME
!     set_orbit_config
!
! DESCRIPTION
!     Setter for 'orbit' configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: real2str, bool2str
use display, only: print_msg

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

! ARGUMENTS
! ---------
logical(k4), intent(in) :: evol_orbit_in, evol_obl_in, evol_ecc_in, evol_lsp_in
real(dp),    intent(in) :: max_change_obl_in, max_change_ecc_in, max_change_lsp_in

! CODE
! ----
evol_orbit = evol_orbit_in
evol_obl = evol_obl_in
evol_ecc = evol_ecc_in
evol_lsp = evol_lsp_in
max_change_obl = max_change_obl_in
max_change_ecc = max_change_ecc_in
max_change_lsp = max_change_lsp_in
call print_msg('evol_orbit     = '//bool2str(evol_orbit))
call print_msg('evol_obl       = '//bool2str(evol_obl))
call print_msg('evol_ecc       = '//bool2str(evol_ecc))
call print_msg('evol_lsp       = '//bool2str(evol_lsp))
call print_msg('max_change_obl = '//real2str(max_change_obl))
call print_msg('max_change_ecc = '//real2str(max_change_ecc))
call print_msg('max_change_lsp = '//real2str(max_change_lsp))

END SUBROUTINE set_orbit_config
!=======================================================================

!=======================================================================
SUBROUTINE init_orbit(obliquity_in,perihelion_in,aphelion_in,date_peri_in,yr_len_sol_in,sol_len_s_in)
!-----------------------------------------------------------------------
! NAME
!     init_orbit
!
! DESCRIPTION
!     Initialize orbital characteristics of the planet.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!     Taken from 'iniorbit' in the Mars PCM.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use maths,   only: pi
use display, only: print_msg
use utility, only: real2str

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: obliquity_in, perihelion_in, aphelion_in, date_peri_in, yr_len_sol_in, sol_len_s_in

! LOCAL VARIABLES
! ---------------
real(dp)    :: zxref, zanom, zz, zx0, zdx
real(dp)    :: lsp ! Solar longitude of the perihelion [rad]
integer(di) :: iter

! CODE
! ----
! Setters
obliquity = obliquity_in
perihelion = perihelion_in
aphelion = aphelion_in
date_peri = date_peri_in
sol_len_s = sol_len_s_in
yr_len_sol = yr_len_sol_in

! Compute eccentricity and semi-major axis
eccentricity = (aphelion - perihelion)/(aphelion + perihelion)
semimaj_axis = 0.5_dp*(aphelion + perihelion)

! Compute mean anomaly zanom
zz = (yr_len_sol - date_peri)/yr_len_sol
zanom = 2._dp*pi*(zz - nint(zz))
zxref = abs(zanom)

! Solve equation zx0 - e*sin(zx0) = zxref for eccentric anomaly zx0 using Newton method
zx0 = zxref + eccentricity*sin(zxref)
do iter = 1,100
    zdx = -(zx0 - eccentricity*sin(zx0) - zxref)/(1._dp - eccentricity*cos(zx0))
    if (abs(zdx) <= 1.e-12_dp) exit
    zx0 = zx0 + zdx
end do

zx0 = zx0 + zdx
if (zanom < 0.) zx0 = -zx0

! Compute solar longitude of the perihelion
lsp = -2._dp*atan(sqrt((1._dp + eccentricity)/(1._dp - eccentricity))*tan(zx0/2._dp))
lsp = modulo(lsp,2._dp*pi)
ls_peri = lsp*180._dp/pi ! Conversion in degree

call print_msg('Obliquity          [deg] = '//real2str(obliquity))
call print_msg('Eccentricity       [-]   = '//real2str(eccentricity))
call print_msg('Ls of perihelion   [deg] = '//real2str(ls_peri))
call print_msg('Perihelion         [Mkm] = '//real2str(perihelion))
call print_msg('Aphelion           [Mkm] = '//real2str(aphelion))
call print_msg('Semi-major axis    [Mkm] = '//real2str(semimaj_axis))
call print_msg('Date of perihelion [sol] = '//real2str(date_peri))
call print_msg('Length of a sol    [s]   = '//real2str(sol_len_s))
call print_msg('Length of a year   [sol] = '//real2str(yr_len_sol))

END SUBROUTINE init_orbit
!=======================================================================

!=======================================================================
SUBROUTINE lsp2datep(lsp,datep)
!-----------------------------------------------------------------------
! NAME
!     lsp2datep
!
! DESCRIPTION
!     Initialize orbital data of the planet.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!     Taken from 'lsp2solp' in the Mars PCM.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use maths, only: pi

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

! ARGUMENTS
! ---------
real(dp), intent(in)  :: lsp   ! Ls at perihelion [deg]
real(dp), intent(out) :: datep ! Date at perihelion [sol]

! LOCAL VARIABLES
! ---------------
real(dp) :: zx0 ! Eccentric anomaly at Ls = 0

! CODE
! ----
! Compute orbit geometrical parameters
zx0 = -2._dp*atan(tan(0.5_dp*lsp*pi/180._dp)*sqrt((1._dp - eccentricity)/(1._dp + eccentricity)))
if (zx0 <= 0._dp) zx0 = zx0 + 2._dp*pi

! Compute sol at perihelion
datep = yr_len_sol*(1._dp - (zx0 - eccentricity*sin(zx0))/(2._dp*pi))

END SUBROUTINE lsp2datep
!=======================================================================

!=======================================================================
SUBROUTINE ini_orbit()
!-----------------------------------------------------------------------
! NAME
!     ini_orbit
!
! DESCRIPTION
!     Initialize the parameters of module 'orbit'.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

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

! LOCAL VARIABLES
! ---------------
integer(di) :: n ! Number of lines in the file
integer(di) :: ierr, funit
logical(k4) :: here

! CODE
! ----
inquire(file = orbitfile_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//orbitfile_name//'"!',1)

! Get the number of lines for Laskar's orbital file
open(newunit = funit,file = orbitfile_name,status = 'old',action = 'read',iostat = ierr)
if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//orbitfile_name//'"!',ierr)
n = 0
do
    read(funit,*,iostat = ierr)
    if (ierr /= 0) exit
    n = n + 1
end do
close(funit)

! Allocation of module arrays
if (.not. allocated(year_lask)) allocate(year_lask(n))
if (.not. allocated(obl_lask)) allocate(obl_lask(n))
if (.not. allocated(ecc_lask)) allocate(ecc_lask(n))
if (.not. allocated(lsp_lask)) allocate(lsp_lask(n))

END SUBROUTINE ini_orbit
!=======================================================================

!=======================================================================
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 read_orbitpm(n_yr_sim)
!-----------------------------------------------------------------------
! NAME
!     read_orbitpm
!
! DESCRIPTION
!     Read Laskar's orbital file and locate index for
!     closest lower year relative to current simulation year.
!
! AUTHORS & DATE
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use evolution, only: pem_ini_date, r_plnt2earth_yr
use stoppage,  only: stop_clean
use display,   only: print_msg
use utility,   only: real2str, int2str

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: n_yr_sim

! LOCAL VARIABLES
! ---------------
integer(di) :: ierr, funit, i, n
logical(k4) :: here
real(dp)    :: current_year

! CODE
! ----
call print_msg('> Reading "'//orbitfile_name//'"')
! Compute the current year
current_year = pem_ini_date + n_yr_sim
call print_msg('Current year = '//real2str(current_year))

inquire(file = orbitfile_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//orbitfile_name//'"!',1)

! Read the Laskars's orbital file
open(newunit = funit,file = orbitfile_name,status = 'old',action = 'read',iostat = ierr)
if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//orbitfile_name//'"!',ierr)
n = size(year_lask,1) ! Number of lines in the file
iyear_lask = 0 ! Line number for closest lower year relative to current simulation year
do i = 1,n
    read(funit,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i)
    year_lask(i) = year_lask(i)*1000._dp/r_plnt2earth_yr ! Conversion from Laskar's kilo Earth years to planetary years
    if (year_lask(i) > current_year) iyear_lask = i + 1
end do
close(funit)
if (iyear_lask == 0 .or. iyear_lask == n + 1) then
    call stop_clean(__FILE__,__LINE__,'the current year could not be found in the file "'//orbitfile_name//'".',1)
else
    call print_msg('The current year in "'//orbitfile_name//'" is at line '//int2str(iyear_lask)//'.')
end if

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

!=======================================================================
SUBROUTINE compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb)
!-----------------------------------------------------------------------
! NAME
!     compute_maxyr_orbit
!
! DESCRIPTION
!     Determine maximum number of years before orbital params change
!     beyond admissible thresholds.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!     Handles Lsp modulo crossings.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use evolution, only: pem_ini_date
use display,   only: print_msg
use utility,   only: real2str

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

! ARGUMENTS
! ---------
real(dp), intent(in)  :: n_yr_sim       ! Current year of the simulation
real(dp), intent(out) :: nmax_yr_runorb ! Maximum number of years of the PEM before orbital parameters change too much

! LOCAL VARIABLES
! ---------------
real(dp)                            :: current_year
real(dp)                            :: max_obl, max_ecc, max_lsp
real(dp)                            :: min_obl, min_ecc, min_lsp
real(dp)                            :: nmax_yr_runobl, nmax_yr_runecc, nmax_yr_runlsp ! Maximum year iteration before reaching an inadmissible value of orbit param
real(dp)                            :: xa, xb, ya, yb
logical(k4)                         :: found_obl, found_ecc, found_lsp
real(dp)                            :: lsp_corr ! Lsp corrected if the "modulo is crossed"
real(dp)                            :: delta
real(dp), dimension(:), allocatable :: lsplask_unwrap
integer(di)                         :: i

! CODE
! ----
if (.not. evol_orbit) then
    nmax_yr_runorb = largest_nb ! Default huge value
    return
end if
call print_msg('> Computing the accepted maximum number of years due to orbital parameters')
current_year = pem_ini_date + n_yr_sim ! We compute the current year

! 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._dp) then
        lsp_corr = lsp_lask(i - 1) - 360._dp
    else if (delta < -180._dp) then
        lsp_corr = lsp_lask(i - 1) + 360._dp
    else
        lsp_corr = lsp_lask(i - 1)
    end if
    lsplask_unwrap(i - 1) = lsplask_unwrap(i) + lsp_corr - lsp_lask(i)
end do

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

! Maximum change accepted for orbital parameters
max_obl = obliquity + max_change_obl
min_obl = obliquity - max_change_obl
call print_msg('Obl. (current|accepted min|accepted max): '//real2str(obliquity)//' | '//real2str(min_obl)//real2str(max_obl))

max_ecc = min(eccentricity + max_change_ecc,1. - minieps) ! ecc < 1.
min_ecc = max(eccentricity - max_change_ecc,0._dp) ! ecc >= 0.
call print_msg('Ecc. (current|accepted min|accepted max): '//real2str(eccentricity)//' | '//real2str(min_ecc)//' | '//real2str(max_ecc))

max_lsp = ls_peri + max_change_lsp
min_lsp = ls_peri - max_change_lsp
call print_msg('Lsp  (current|accepted min|accepted max): '//real2str(ls_peri)//' | '//real2str(min_lsp)//' | '//real2str(max_lsp))

! Initialization
nmax_yr_runobl = largest_nb ! Default huge value
nmax_yr_runecc = largest_nb ! Default huge value
nmax_yr_runlsp =largest_nb ! Default huge value
found_obl = .false.
found_ecc = .false.
found_lsp = .false.
if (.not. evol_obl) found_obl = .true.
if (.not. evol_ecc) found_ecc = .true.
if (.not. evol_lsp) found_lsp = .true.

! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
do i = iyear_lask,2,-1
    xa = year_lask(i)
    xb = year_lask(i - 1)

    ! Obliquity
    if (evol_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
            nmax_yr_runobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - current_year
            found_obl = .true.
        else if (max_obl < yb) then ! The maximum accepted is overtaken
            nmax_yr_runobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - current_year
            found_obl = .true.
        end if
    end if

    ! Eccentricity
    if (evol_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
            nmax_yr_runecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - current_year
            found_ecc = .true.
        else if (max_ecc < yb) then ! The maximum accepted is overtaken
            nmax_yr_runecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - current_year
            found_ecc = .true.
        end if
    end if

    ! Lsp
    if (evol_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
            nmax_yr_runlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - current_year
            found_lsp = .true.
        else if (max_lsp < yb) then ! The maximum accepted is overtaken
            nmax_yr_runlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - current_year
            found_lsp = .true.
        end if
    end if

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

deallocate(lsplask_unwrap)

if (.not. found_obl) call print_msg('The maximum reachable year for obl. could not be found in "'//orbitfile_name//'".')
if (.not. found_ecc) call print_msg('The maximum reachable year for ecc. could not be found in "'//orbitfile_name//'".')
if (.not. found_lsp) call print_msg('The maximum reachable year for Lsp could not be found in "'//orbitfile_name//'".')

call print_msg('Number of years accepted for obl. = '//real2str(nmax_yr_runobl))
call print_msg('Number of years accepted for ecc. = '//real2str(nmax_yr_runecc))
call print_msg('Number of years accepted for Lsp  = '//real2str(nmax_yr_runlsp))

nmax_yr_runorb = min(nmax_yr_runobl,nmax_yr_runecc,nmax_yr_runlsp)
if (nmax_yr_runorb < 1._dp) nmax_yr_runorb = 1._dp ! nmax_yr_runorb should be at least equal to 1

call print_msg('So the max. number of years for the orbital criterion = '//real2str(nmax_yr_runorb))

END SUBROUTINE compute_maxyr_orbit
!=======================================================================

!=======================================================================
SUBROUTINE update_orbit(n_yr_sim,n_yr_run)
!-----------------------------------------------------------------------
! NAME
!     update_orbit
!
! DESCRIPTION
!     Update orbital parameters from Laskar values at a new year.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

use evolution,        only: pem_ini_date
use call_dayperi_mod, only: call_dayperi
use stoppage,         only: stop_clean
use display,          only: print_msg
use utility,          only: real2str

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: n_yr_sim ! Number of simulated years
real(dp), intent(in) :: n_yr_run ! Number of iterations of the PEM

! LOCAL VARIABLES
! ---------------
integer(di) :: i
real(dp)    :: current_year, old_year, obl_old, ecc_old, lsp_old
real(dp)    :: xa, xb, ya, yb ! Variables for interpolation
logical(k4) :: found_year

! CODE
! ----
call print_msg('> Updating the orbital parameters for the new year')
current_year = pem_ini_date + n_yr_sim ! We compute the current year
old_year = current_year - n_yr_run
obl_old = obliquity
ecc_old = eccentricity
lsp_old = ls_peri

found_year = .false.
if (current_year < 0._dp) then
    do i = iyear_lask,2,-1
        xa = year_lask(i)
        xb = year_lask(i - 1)
        if (xa <= current_year .and. current_year < xb) then
            ! Obliquity
            if (evol_obl) then
                ya = obl_lask(i)
                yb = obl_lask(i - 1)
                obliquity = (current_year - xa)*(yb - ya)/(xb - xa) + ya
            end if

            ! Eccentricity
            if (evol_ecc) then
                ya = ecc_lask(i)
                yb = ecc_lask(i - 1)
                eccentricity = (current_year - xa)*(yb - ya)/(xb - xa) + ya
            end if

            ! Lsp
            if (evol_lsp) then
                ya = lsp_lask(i)
                yb = lsp_lask(i - 1)
                if (abs(yb - ya) > 300._dp) then ! If modulo is "crossed" through the interval
                    if (ya < yb) then ! Lsp should be decreasing
                        ya = ya + 360._dp
                    else ! Lsp should be increasing
                        yb = yb + 360._dp
                    end if
                end if
                ls_peri = modulo((current_year - xa)*(yb - ya)/(xb - xa) + ya,360._dp)
            end if
            found_year = .true.
            exit ! The loop is left as soon as the right interval is found
        end if
    end do
else if (current_year >= 0._dp .and. current_year < 100._dp) then ! For present orbital characteristics
    if (evol_obl) obliquity = obl_lask(1)
    if (evol_ecc) eccentricity = ecc_lask(1)
    if (evol_lsp) ls_peri = lsp_lask(1)
    found_year = .true.
end if

if (.not. found_year) call stop_clean(__FILE__,__LINE__,'the new year could not be found in the file "'//orbitfile_name//'".',1)

call lsp2datep(ls_peri,date_peri)
perihelion = semimaj_axis*(1._dp - eccentricity)
aphelion = semimaj_axis*(1._dp + eccentricity)

call print_msg('Year (ini|now) = '//real2str(old_year)//' | '//real2str(current_year))
call print_msg('Obl. (ini|now) = '//real2str(obl_old)//' | '//real2str(obliquity))
call print_msg('Ecc. (ini|now) = '//real2str(ecc_old)//' | '//real2str(eccentricity))
call print_msg('Lsp  (ini|now) = '//real2str(lsp_old)//' | '//real2str(ls_peri))

END SUBROUTINE update_orbit
!=======================================================================

END MODULE orbit