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

MODULE orbit

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

implicit none

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)

use evolution, only: convert_years

implicit none

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

integer :: n ! number of lines in Laskar files
integer :: ierr, i

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

implicit none

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)

#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

implicit none

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

!--------------------------------------------------------
! Output Variables
!--------------------------------------------------------
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
integer                         :: i                                          ! Loop variable
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

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

!***********************************************************************
! Purpose: Recompute orbit parameters based on values by Laskar et al.
SUBROUTINE set_new_orbitpm(i_myear,i_myear_leg)

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

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 :: i          ! Loop variable
real    :: halfaxe        ! Million of km
real    :: unitastr
real    :: xa, xb, ya, yb ! Variables for interpolation
logical :: found_year     ! Flag variable

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