source: trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90 @ 3985

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,n_myear_leg)

#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 time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask

implicit none

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

!--------------------------------------------------------
! Output Variables
!--------------------------------------------------------
real, intent(out) :: n_myear_leg ! 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                            :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
integer                         :: ilask                                          ! 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 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',status = 'old',action = 'read')
last_ilask = 0
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)
if (last_ilask == 0 .or. last_ilask == size(yearlask,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:', last_ilask
endif

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

! Correct the current Lsp according to the unwrapping
delta = Lsp - lsplask_unwrap(last_ilask)
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
max_obl_iter = 1000000.
max_ecc_iter = 1000000.
max_lsp_iter = 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 ilask = last_ilask,2,-1
    xa = yearlask(ilask)
    xb = yearlask(ilask - 1)

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

    ! Eccentricity
    if (var_ecc .and. .not. found_ecc) then
        ya = ecclask(ilask)
        yb = ecclask(ilask - 1)
        if (yb < min_ecc) then ! The minimum accepted is overtaken
            max_ecc_iter = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_ecc = .true.
        else if (max_ecc < yb) then ! The maximum accepted is overtaken
            max_ecc_iter = (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(ilask)
        yb = lsplask_unwrap(ilask - 1)
        if (yb < min_lsp) then ! The minimum accepted is overtaken
            max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
            found_lsp = .true.
        else if (max_lsp < yb) then ! The maximum accepted is overtaken
            max_lsp_iter = (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 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

n_myear_leg = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
if (n_myear_leg < 1.) n_myear_leg = 1. ! n_myear_leg should be at least equal to 1
write(*,*) 'So the max. number of iterations for the orbital criteria =', n_myear_leg

END SUBROUTINE orbit_param_criterion

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

END MODULE orbit_param_criterion_mod