| 1 | MODULE orbit_param_criterion_mod |
|---|
| 2 | |
|---|
| 3 | !======================================================================= |
|---|
| 4 | ! |
|---|
| 5 | ! Purpose: Compute the maximum number of iteration for PEM based on the |
|---|
| 6 | ! stopping criterion given by the orbital parameters |
|---|
| 7 | ! |
|---|
| 8 | ! Author: RV, JBC |
|---|
| 9 | !======================================================================= |
|---|
| 10 | |
|---|
| 11 | implicit none |
|---|
| 12 | |
|---|
| 13 | !======================================================================= |
|---|
| 14 | contains |
|---|
| 15 | !======================================================================= |
|---|
| 16 | |
|---|
| 17 | SUBROUTINE orbit_param_criterion(i_myear,n_myear_leg) |
|---|
| 18 | |
|---|
| 19 | #ifdef CPP_IOIPSL |
|---|
| 20 | use IOIPSL, only: getin |
|---|
| 21 | #else |
|---|
| 22 | ! if not using IOIPSL, we still need to use (a local version of) getin |
|---|
| 23 | use ioipsl_getincom, only: getin |
|---|
| 24 | #endif |
|---|
| 25 | #ifndef CPP_STD |
|---|
| 26 | use planete_h, only: e_elips, obliquit, lsperi |
|---|
| 27 | use comcstfi_h, only: pi |
|---|
| 28 | #else |
|---|
| 29 | use planete_mod, only: e_elips, obliquit, lsperi |
|---|
| 30 | use comcstfi_mod, only: pi |
|---|
| 31 | #endif |
|---|
| 32 | use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years |
|---|
| 33 | use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask |
|---|
| 34 | |
|---|
| 35 | implicit none |
|---|
| 36 | |
|---|
| 37 | !-------------------------------------------------------- |
|---|
| 38 | ! Input Variables |
|---|
| 39 | !-------------------------------------------------------- |
|---|
| 40 | real, intent(in) :: i_myear ! Martian year of the simulation |
|---|
| 41 | |
|---|
| 42 | !-------------------------------------------------------- |
|---|
| 43 | ! Output Variables |
|---|
| 44 | !-------------------------------------------------------- |
|---|
| 45 | real, intent(out) :: n_myear_leg ! Maximum number of iteration of the PEM before orb changes too much |
|---|
| 46 | |
|---|
| 47 | !-------------------------------------------------------- |
|---|
| 48 | ! Local variables |
|---|
| 49 | !-------------------------------------------------------- |
|---|
| 50 | real :: Year ! Year of the simulation |
|---|
| 51 | real :: Lsp ! Ls perihelion |
|---|
| 52 | real :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change |
|---|
| 53 | real :: max_obl, max_ecc, max_lsp ! Maximum value of orbit param given the admissible change |
|---|
| 54 | real :: min_obl, min_ecc, min_lsp ! Minimum value of orbit param given the admissible change |
|---|
| 55 | real :: max_obl_iter, max_ecc_iter, max_lsp_iter ! Maximum year iteration before reaching an inadmissible value of orbit param |
|---|
| 56 | integer :: ilask ! Loop variable |
|---|
| 57 | real :: xa, xb, ya, yb ! Variables for interpolation |
|---|
| 58 | logical :: found_obl, found_ecc, found_lsp ! Flag variables for orbital parameters |
|---|
| 59 | real :: lsp_corr ! Lsp corrected if the "modulo is crossed" |
|---|
| 60 | real :: delta ! Intermediate variable |
|---|
| 61 | real, dimension(:), allocatable :: lsplask_unwrap ! Unwrapped sequence of Lsp from Laskar's file |
|---|
| 62 | |
|---|
| 63 | ! ********************************************************************** |
|---|
| 64 | ! 0. Initializations |
|---|
| 65 | ! ********************************************************************** |
|---|
| 66 | Year = year_bp_ini + i_myear ! We compute the current year |
|---|
| 67 | Lsp = lsperi*180./pi ! We convert in degree |
|---|
| 68 | |
|---|
| 69 | call ini_lask_param_mod() ! Allocation of variables |
|---|
| 70 | |
|---|
| 71 | write(*,*) 'orbit_param_criterion, Current year =', Year |
|---|
| 72 | write(*,*) 'orbit_param_criterion, Current obl. =', obliquit |
|---|
| 73 | write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips |
|---|
| 74 | write(*,*) 'orbit_param_criterion, Current Lsp =', Lsp |
|---|
| 75 | |
|---|
| 76 | ! We read the file |
|---|
| 77 | open(73,file = 'obl_ecc_lsp.asc') |
|---|
| 78 | last_ilask = 0 |
|---|
| 79 | do ilask = 1,size(yearlask,1) |
|---|
| 80 | read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask) |
|---|
| 81 | yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years |
|---|
| 82 | if (yearlask(ilask) > Year) last_ilask = ilask + 1 |
|---|
| 83 | enddo |
|---|
| 84 | close(73) |
|---|
| 85 | if (last_ilask == 0 .or. last_ilask == size(yearlask,1) + 1) then |
|---|
| 86 | error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 87 | else |
|---|
| 88 | write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', last_ilask |
|---|
| 89 | endif |
|---|
| 90 | |
|---|
| 91 | ! Unwrap the Lsp changes in Laskar's file |
|---|
| 92 | allocate(lsplask_unwrap(last_ilask)) |
|---|
| 93 | lsplask_unwrap(last_ilask) = lsplask(last_ilask) |
|---|
| 94 | do ilask = last_ilask,2,-1 |
|---|
| 95 | delta = lsplask(ilask - 1) - lsplask(ilask) |
|---|
| 96 | if (delta > 180.) then |
|---|
| 97 | lsp_corr = lsplask(ilask - 1) - 360. |
|---|
| 98 | else if (delta < -180.) then |
|---|
| 99 | lsp_corr = lsplask(ilask - 1) + 360. |
|---|
| 100 | else |
|---|
| 101 | lsp_corr = lsplask(ilask - 1) |
|---|
| 102 | endif |
|---|
| 103 | lsplask_unwrap(ilask - 1) = lsplask_unwrap(ilask) + lsp_corr - lsplask(ilask) |
|---|
| 104 | enddo |
|---|
| 105 | |
|---|
| 106 | ! Correct the current Lsp according to the unwrapping |
|---|
| 107 | delta = Lsp - lsplask_unwrap(last_ilask) |
|---|
| 108 | if (delta > 180.) then |
|---|
| 109 | Lsp = Lsp - 360. |
|---|
| 110 | else if (delta < -180.) then |
|---|
| 111 | Lsp = Lsp + 360. |
|---|
| 112 | endif |
|---|
| 113 | |
|---|
| 114 | ! Constant max change case |
|---|
| 115 | max_change_obl = 1. |
|---|
| 116 | max_change_ecc = 5.e-3 |
|---|
| 117 | max_change_lsp = 20. |
|---|
| 118 | |
|---|
| 119 | call getin('max_change_obl',max_change_obl) |
|---|
| 120 | max_obl = obliquit + max_change_obl |
|---|
| 121 | min_obl = obliquit - max_change_obl |
|---|
| 122 | write(*,*) 'Maximum obliquity accepted =', max_obl |
|---|
| 123 | write(*,*) 'Minimum obliquity accepted =', min_obl |
|---|
| 124 | |
|---|
| 125 | call getin('max_change_ecc',max_change_ecc) |
|---|
| 126 | max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ecc < 1. |
|---|
| 127 | min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 0. |
|---|
| 128 | write(*,*) 'Maximum eccentricity accepted =', max_ecc |
|---|
| 129 | write(*,*) 'Minimum eccentricity accepted =', min_ecc |
|---|
| 130 | |
|---|
| 131 | call getin('max_change_lsp',max_change_lsp) |
|---|
| 132 | max_lsp = Lsp + max_change_lsp |
|---|
| 133 | min_lsp = Lsp - max_change_lsp |
|---|
| 134 | write(*,*) 'Maximum Lsp accepted =', max_lsp |
|---|
| 135 | write(*,*) 'Minimum Lsp accepted =', min_lsp |
|---|
| 136 | ! End Constant max change case |
|---|
| 137 | |
|---|
| 138 | ! If we do not want some orb parameter to change, they should not be a stopping criterion, |
|---|
| 139 | ! So the number of iteration corresponding is set to maximum |
|---|
| 140 | max_obl_iter = 1000000. |
|---|
| 141 | max_ecc_iter = 1000000. |
|---|
| 142 | max_lsp_iter = 1000000. |
|---|
| 143 | |
|---|
| 144 | ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters |
|---|
| 145 | found_obl = .false. |
|---|
| 146 | found_ecc = .false. |
|---|
| 147 | found_lsp = .false. |
|---|
| 148 | if (.not. var_obl) found_obl = .true. |
|---|
| 149 | if (.not. var_ecc) found_ecc = .true. |
|---|
| 150 | if (.not. var_lsp) found_lsp = .true. |
|---|
| 151 | |
|---|
| 152 | do ilask = last_ilask,2,-1 |
|---|
| 153 | xa = yearlask(ilask) |
|---|
| 154 | xb = yearlask(ilask - 1) |
|---|
| 155 | |
|---|
| 156 | ! Obliquity |
|---|
| 157 | if (var_obl .and. .not. found_obl) then |
|---|
| 158 | ya = obllask(ilask) |
|---|
| 159 | yb = obllask(ilask - 1) |
|---|
| 160 | if (yb < min_obl) then ! The minimum accepted is overtaken |
|---|
| 161 | max_obl_iter = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 162 | found_obl = .true. |
|---|
| 163 | else if (max_obl < yb) then ! The maximum accepted is overtaken |
|---|
| 164 | max_obl_iter = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 165 | found_obl = .true. |
|---|
| 166 | endif |
|---|
| 167 | endif |
|---|
| 168 | |
|---|
| 169 | ! Eccentricity |
|---|
| 170 | if (var_ecc .and. .not. found_ecc) then |
|---|
| 171 | ya = ecclask(ilask) |
|---|
| 172 | yb = ecclask(ilask - 1) |
|---|
| 173 | if (yb < min_ecc) then ! The minimum accepted is overtaken |
|---|
| 174 | max_ecc_iter = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 175 | found_ecc = .true. |
|---|
| 176 | else if (max_ecc < yb) then ! The maximum accepted is overtaken |
|---|
| 177 | max_ecc_iter = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 178 | found_ecc = .true. |
|---|
| 179 | endif |
|---|
| 180 | endif |
|---|
| 181 | |
|---|
| 182 | ! Lsp |
|---|
| 183 | if (var_lsp .and. .not. found_lsp) then |
|---|
| 184 | ya = lsplask_unwrap(ilask) |
|---|
| 185 | yb = lsplask_unwrap(ilask - 1) |
|---|
| 186 | if (yb < min_lsp) then ! The minimum accepted is overtaken |
|---|
| 187 | max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 188 | found_lsp = .true. |
|---|
| 189 | else if (max_lsp < yb) then ! The maximum accepted is overtaken |
|---|
| 190 | max_lsp_iter = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 191 | found_lsp = .true. |
|---|
| 192 | endif |
|---|
| 193 | endif |
|---|
| 194 | |
|---|
| 195 | if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found |
|---|
| 196 | enddo |
|---|
| 197 | |
|---|
| 198 | deallocate(lsplask_unwrap) |
|---|
| 199 | |
|---|
| 200 | if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 201 | if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 202 | if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 203 | |
|---|
| 204 | write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter |
|---|
| 205 | write(*,*) 'Max. number of iterations for the ecc. parameter =', max_ecc_iter |
|---|
| 206 | write(*,*) 'Max. number of iterations for the Lsp parameter =', max_lsp_iter |
|---|
| 207 | |
|---|
| 208 | n_myear_leg = min(max_obl_iter,max_ecc_iter,max_lsp_iter) |
|---|
| 209 | if (n_myear_leg < 1.) n_myear_leg = 1. ! n_myear_leg should be at least equal to 1 |
|---|
| 210 | write(*,*) 'So the max. number of iterations for the orbital criteria =', n_myear_leg |
|---|
| 211 | |
|---|
| 212 | END SUBROUTINE orbit_param_criterion |
|---|
| 213 | |
|---|
| 214 | !*********************************************************************** |
|---|
| 215 | |
|---|
| 216 | END MODULE orbit_param_criterion_mod |
|---|
| 217 | |
|---|
