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 #ifndef CPP_STD use planete_h, only: e_elips, obliquit, lsperi use comcstfi_h, only: pi #else use planete_mod, only: e_elips, obliquit, lsperi use comcstfi_mod, only: pi #endif 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') 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