source: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 @ 3860

MODULE recomp_orb_param_mod
!=======================================================================
!
! Purpose: Recompute orbit parameters based on values by Laskar et al.,
!
! Author: RV, JBC
!=======================================================================

implicit none

!=======================================================================
contains
!=======================================================================

SUBROUTINE recomp_orb_param(i_myear,i_myear_leg)

use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp
use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
#ifndef CPP_STD
    use comcstfi_h,       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
#else
    use comcstfi_mod,     only: pi
    use planete_mod,      only: e_elips, obliquit, lsperi, periastr, apoastr, p_elips, peri_day, year_day
#endif

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

! **********************************************************************
! 0. Initializations
! **********************************************************************
#ifdef CPP_STD
    real :: aphelie, periheli

    aphelie = apoastr
    periheli = periastr
#endif

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

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

found_year = .false.
if (Year < 0.) then
    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 = obllask(ilask)
                yb = obllask(ilask - 1)
                obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
            endif

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

            ! Lsp
            if (var_lsp) then
                ya = lsplask(ilask)
                yb = lsplask(ilask - 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(*,*) 'recomp_orb_param, New obl. =', obliquit
write(*,*) 'recomp_orb_param, New ecc. =', e_elips
write(*,*) 'recomp_orb_param, 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

#ifndef CPP_STD
    call call_dayperi(Lsp,e_elips,peri_day,year_day)
#endif

END SUBROUTINE recomp_orb_param

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

END MODULE recomp_orb_param_mod