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