| 1 | MODULE orbit |
|---|
| 2 | !----------------------------------------------------------------------- |
|---|
| 3 | ! NAME |
|---|
| 4 | ! orbit |
|---|
| 5 | ! DESCRIPTION |
|---|
| 6 | ! Compute orbital-parameter-based stopping criteria and update |
|---|
| 7 | ! planetary orbit parameters. |
|---|
| 8 | ! AUTHORS & DATE |
|---|
| 9 | ! R. Vandemeulebrouck |
|---|
| 10 | ! JB Clement, 2023-2025 |
|---|
| 11 | ! NOTES |
|---|
| 12 | ! |
|---|
| 13 | !----------------------------------------------------------------------- |
|---|
| 14 | |
|---|
| 15 | ! DECLARATION |
|---|
| 16 | ! ----------- |
|---|
| 17 | implicit none |
|---|
| 18 | |
|---|
| 19 | ! MODULE VARIABLES |
|---|
| 20 | ! ---------------- |
|---|
| 21 | real, dimension(:), allocatable :: year_lask ! year before present from Laskar et al. |
|---|
| 22 | real, dimension(:), allocatable :: obl_lask ! obliquity [deg] |
|---|
| 23 | real, dimension(:), allocatable :: ecc_lask ! eccentricity |
|---|
| 24 | real, dimension(:), allocatable :: lsp_lask ! ls perihelie [deg] |
|---|
| 25 | integer :: iyear_lask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year |
|---|
| 26 | |
|---|
| 27 | contains |
|---|
| 28 | !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|---|
| 29 | |
|---|
| 30 | !======================================================================= |
|---|
| 31 | SUBROUTINE read_orbitpm(Year) |
|---|
| 32 | !----------------------------------------------------------------------- |
|---|
| 33 | ! NAME |
|---|
| 34 | ! read_orbitpm |
|---|
| 35 | ! |
|---|
| 36 | ! DESCRIPTION |
|---|
| 37 | ! Read Laskar orbital file and prepare arrays; locate index for |
|---|
| 38 | ! closest lower year relative to current simulation year. |
|---|
| 39 | ! |
|---|
| 40 | ! AUTHORS & DATE |
|---|
| 41 | ! R. Vandemeulebrouck |
|---|
| 42 | ! JB Clement, 2023-2025 |
|---|
| 43 | ! |
|---|
| 44 | ! NOTES |
|---|
| 45 | ! Converts kilo Earth years to Martian years using 'convert_years'. |
|---|
| 46 | !----------------------------------------------------------------------- |
|---|
| 47 | |
|---|
| 48 | ! DEPENDENCIES |
|---|
| 49 | ! ------------ |
|---|
| 50 | use evolution, only: convert_years |
|---|
| 51 | |
|---|
| 52 | ! DECLARATION |
|---|
| 53 | ! ----------- |
|---|
| 54 | implicit none |
|---|
| 55 | |
|---|
| 56 | ! ARGUMENTS |
|---|
| 57 | ! --------- |
|---|
| 58 | real, intent(in) :: Year ! Martian year of the simulation |
|---|
| 59 | |
|---|
| 60 | ! LOCAL VARIABLES |
|---|
| 61 | ! --------------- |
|---|
| 62 | integer :: n ! Number of lines in Laskar files |
|---|
| 63 | integer :: ierr, i |
|---|
| 64 | |
|---|
| 65 | ! CODE |
|---|
| 66 | ! ---- |
|---|
| 67 | n = 0 |
|---|
| 68 | open(1,file = 'obl_ecc_lsp.asc',status = 'old',action = 'read') |
|---|
| 69 | do |
|---|
| 70 | read(1,*,iostat = ierr) |
|---|
| 71 | if (ierr /= 0) exit |
|---|
| 72 | n = n + 1 |
|---|
| 73 | enddo |
|---|
| 74 | close(1) |
|---|
| 75 | allocate(year_lask(n),obl_lask(n),ecc_lask(n),lsp_lask(n)) |
|---|
| 76 | |
|---|
| 77 | iyear_lask = 0 |
|---|
| 78 | do i = 1,size(year_lask,1) |
|---|
| 79 | read(1,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i) |
|---|
| 80 | year_lask(i) = year_lask(i)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years |
|---|
| 81 | if (year_lask(i) > Year) iyear_lask = i + 1 |
|---|
| 82 | enddo |
|---|
| 83 | close(1) |
|---|
| 84 | if (iyear_lask == 0 .or. iyear_lask == size(year_lask,1) + 1) then |
|---|
| 85 | error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 86 | else |
|---|
| 87 | write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', iyear_lask |
|---|
| 88 | endif |
|---|
| 89 | |
|---|
| 90 | END SUBROUTINE read_orbitpm |
|---|
| 91 | !======================================================================= |
|---|
| 92 | |
|---|
| 93 | !======================================================================= |
|---|
| 94 | SUBROUTINE end_orbit |
|---|
| 95 | !----------------------------------------------------------------------- |
|---|
| 96 | ! NAME |
|---|
| 97 | ! end_orbit |
|---|
| 98 | ! |
|---|
| 99 | ! DESCRIPTION |
|---|
| 100 | ! Deallocate orbital data arrays. |
|---|
| 101 | ! |
|---|
| 102 | ! AUTHORS & DATE |
|---|
| 103 | ! R. Vandemeulebrouck |
|---|
| 104 | ! JB Clement, 2023-2025 |
|---|
| 105 | ! |
|---|
| 106 | ! NOTES |
|---|
| 107 | ! |
|---|
| 108 | !----------------------------------------------------------------------- |
|---|
| 109 | |
|---|
| 110 | ! DECLARATION |
|---|
| 111 | ! ----------- |
|---|
| 112 | implicit none |
|---|
| 113 | |
|---|
| 114 | ! CODE |
|---|
| 115 | ! ---- |
|---|
| 116 | if (allocated(year_lask)) deallocate(year_lask) |
|---|
| 117 | if (allocated(obl_lask)) deallocate(obl_lask) |
|---|
| 118 | if (allocated(ecc_lask)) deallocate(ecc_lask) |
|---|
| 119 | if (allocated(lsp_lask)) deallocate(lsp_lask) |
|---|
| 120 | |
|---|
| 121 | END SUBROUTINE end_orbit |
|---|
| 122 | !======================================================================= |
|---|
| 123 | |
|---|
| 124 | !======================================================================= |
|---|
| 125 | SUBROUTINE orbit_param_criterion(i_myear,nyears_maxorb) |
|---|
| 126 | !----------------------------------------------------------------------- |
|---|
| 127 | ! NAME |
|---|
| 128 | ! orbit_param_criterion |
|---|
| 129 | ! |
|---|
| 130 | ! DESCRIPTION |
|---|
| 131 | ! Determine maximum PEM iterations before orbital params change |
|---|
| 132 | ! beyond admissible thresholds. |
|---|
| 133 | ! |
|---|
| 134 | ! AUTHORS & DATE |
|---|
| 135 | ! R. Vandemeulebrouck |
|---|
| 136 | ! JB Clement, 2023-2025 |
|---|
| 137 | ! |
|---|
| 138 | ! NOTES |
|---|
| 139 | ! |
|---|
| 140 | !----------------------------------------------------------------------- |
|---|
| 141 | #ifdef CPP_IOIPSL |
|---|
| 142 | use IOIPSL, only: getin |
|---|
| 143 | #else |
|---|
| 144 | ! if not using IOIPSL, we still need to use (a local version of) getin |
|---|
| 145 | use ioipsl_getincom, only: getin |
|---|
| 146 | #endif |
|---|
| 147 | use planete_h, only: e_elips, obliquit, lsperi |
|---|
| 148 | use phys_constants, only: pi |
|---|
| 149 | use evolution, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years |
|---|
| 150 | |
|---|
| 151 | ! DECLARATION |
|---|
| 152 | ! ----------- |
|---|
| 153 | implicit none |
|---|
| 154 | |
|---|
| 155 | ! ARGUMENTS |
|---|
| 156 | ! --------- |
|---|
| 157 | real, intent(in) :: i_myear ! Martian year of the simulation |
|---|
| 158 | real, intent(out) :: nyears_maxorb ! Maximum number of iteration of the PEM before orb changes too much |
|---|
| 159 | |
|---|
| 160 | ! LOCAL VARIABLES |
|---|
| 161 | ! --------------- |
|---|
| 162 | real :: Year ! Year of the simulation |
|---|
| 163 | real :: Lsp ! Ls perihelion |
|---|
| 164 | real :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change |
|---|
| 165 | real :: max_obl, max_ecc, max_lsp ! Maximum value of orbit param given the admissible change |
|---|
| 166 | real :: min_obl, min_ecc, min_lsp ! Minimum value of orbit param given the admissible change |
|---|
| 167 | real :: nyears_maxobl, nyears_maxecc, nyears_maxlsp ! Maximum year iteration before reaching an inadmissible value of orbit param |
|---|
| 168 | real :: xa, xb, ya, yb ! Variables for interpolation |
|---|
| 169 | logical :: found_obl, found_ecc, found_lsp ! Flag variables for orbital parameters |
|---|
| 170 | real :: lsp_corr ! Lsp corrected if the "modulo is crossed" |
|---|
| 171 | real :: delta ! Intermediate variable |
|---|
| 172 | real, dimension(:), allocatable :: lsplask_unwrap ! Unwrapped sequence of Lsp from Laskar's file |
|---|
| 173 | integer :: i |
|---|
| 174 | |
|---|
| 175 | ! CODE |
|---|
| 176 | ! ---- |
|---|
| 177 | ! ********************************************************************** |
|---|
| 178 | ! 0. Initializations |
|---|
| 179 | ! ********************************************************************** |
|---|
| 180 | Year = year_bp_ini + i_myear ! We compute the current year |
|---|
| 181 | Lsp = lsperi*180./pi ! We convert in degree |
|---|
| 182 | |
|---|
| 183 | call read_orbitpm(Year) ! Allocation of variables |
|---|
| 184 | |
|---|
| 185 | write(*,*) 'orbit_param_criterion, Current year =', Year |
|---|
| 186 | write(*,*) 'orbit_param_criterion, Current obl. =', obliquit |
|---|
| 187 | write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips |
|---|
| 188 | write(*,*) 'orbit_param_criterion, Current Lsp =', Lsp |
|---|
| 189 | |
|---|
| 190 | ! Unwrap the Lsp changes in Laskar's file |
|---|
| 191 | allocate(lsplask_unwrap(iyear_lask)) |
|---|
| 192 | lsplask_unwrap(iyear_lask) = lsp_lask(iyear_lask) |
|---|
| 193 | do i = iyear_lask,2,-1 |
|---|
| 194 | delta = lsp_lask(i - 1) - lsp_lask(i) |
|---|
| 195 | if (delta > 180.) then |
|---|
| 196 | lsp_corr = lsp_lask(i - 1) - 360. |
|---|
| 197 | else if (delta < -180.) then |
|---|
| 198 | lsp_corr = lsp_lask(i - 1) + 360. |
|---|
| 199 | else |
|---|
| 200 | lsp_corr = lsp_lask(i - 1) |
|---|
| 201 | endif |
|---|
| 202 | lsplask_unwrap(i - 1) = lsplask_unwrap(i) + lsp_corr - lsp_lask(i) |
|---|
| 203 | enddo |
|---|
| 204 | |
|---|
| 205 | ! Correct the current Lsp according to the unwrapping |
|---|
| 206 | delta = Lsp - lsplask_unwrap(iyear_lask) |
|---|
| 207 | if (delta > 180.) then |
|---|
| 208 | Lsp = Lsp - 360. |
|---|
| 209 | else if (delta < -180.) then |
|---|
| 210 | Lsp = Lsp + 360. |
|---|
| 211 | endif |
|---|
| 212 | |
|---|
| 213 | ! Constant max change case |
|---|
| 214 | max_change_obl = 1. |
|---|
| 215 | max_change_ecc = 5.e-3 |
|---|
| 216 | max_change_lsp = 20. |
|---|
| 217 | |
|---|
| 218 | call getin('max_change_obl',max_change_obl) |
|---|
| 219 | max_obl = obliquit + max_change_obl |
|---|
| 220 | min_obl = obliquit - max_change_obl |
|---|
| 221 | write(*,*) 'Maximum obliquity accepted =', max_obl |
|---|
| 222 | write(*,*) 'Minimum obliquity accepted =', min_obl |
|---|
| 223 | |
|---|
| 224 | call getin('max_change_ecc',max_change_ecc) |
|---|
| 225 | max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ecc < 1. |
|---|
| 226 | min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 0. |
|---|
| 227 | write(*,*) 'Maximum eccentricity accepted =', max_ecc |
|---|
| 228 | write(*,*) 'Minimum eccentricity accepted =', min_ecc |
|---|
| 229 | |
|---|
| 230 | call getin('max_change_lsp',max_change_lsp) |
|---|
| 231 | max_lsp = Lsp + max_change_lsp |
|---|
| 232 | min_lsp = Lsp - max_change_lsp |
|---|
| 233 | write(*,*) 'Maximum Lsp accepted =', max_lsp |
|---|
| 234 | write(*,*) 'Minimum Lsp accepted =', min_lsp |
|---|
| 235 | ! End Constant max change case |
|---|
| 236 | |
|---|
| 237 | ! If we do not want some orb parameter to change, they should not be a stopping criterion, |
|---|
| 238 | ! So the number of iteration corresponding is set to maximum |
|---|
| 239 | nyears_maxobl = 1000000. |
|---|
| 240 | nyears_maxecc = 1000000. |
|---|
| 241 | nyears_maxlsp = 1000000. |
|---|
| 242 | |
|---|
| 243 | ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters |
|---|
| 244 | found_obl = .false. |
|---|
| 245 | found_ecc = .false. |
|---|
| 246 | found_lsp = .false. |
|---|
| 247 | if (.not. var_obl) found_obl = .true. |
|---|
| 248 | if (.not. var_ecc) found_ecc = .true. |
|---|
| 249 | if (.not. var_lsp) found_lsp = .true. |
|---|
| 250 | |
|---|
| 251 | do i = iyear_lask,2,-1 |
|---|
| 252 | xa = year_lask(i) |
|---|
| 253 | xb = year_lask(i - 1) |
|---|
| 254 | |
|---|
| 255 | ! Obliquity |
|---|
| 256 | if (var_obl .and. .not. found_obl) then |
|---|
| 257 | ya = obl_lask(i) |
|---|
| 258 | yb = obl_lask(i - 1) |
|---|
| 259 | if (yb < min_obl) then ! The minimum accepted is overtaken |
|---|
| 260 | nyears_maxobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 261 | found_obl = .true. |
|---|
| 262 | else if (max_obl < yb) then ! The maximum accepted is overtaken |
|---|
| 263 | nyears_maxobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 264 | found_obl = .true. |
|---|
| 265 | endif |
|---|
| 266 | endif |
|---|
| 267 | |
|---|
| 268 | ! Eccentricity |
|---|
| 269 | if (var_ecc .and. .not. found_ecc) then |
|---|
| 270 | ya = ecc_lask(i) |
|---|
| 271 | yb = ecc_lask(i - 1) |
|---|
| 272 | if (yb < min_ecc) then ! The minimum accepted is overtaken |
|---|
| 273 | nyears_maxecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 274 | found_ecc = .true. |
|---|
| 275 | else if (max_ecc < yb) then ! The maximum accepted is overtaken |
|---|
| 276 | nyears_maxecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 277 | found_ecc = .true. |
|---|
| 278 | endif |
|---|
| 279 | endif |
|---|
| 280 | |
|---|
| 281 | ! Lsp |
|---|
| 282 | if (var_lsp .and. .not. found_lsp) then |
|---|
| 283 | ya = lsplask_unwrap(i) |
|---|
| 284 | yb = lsplask_unwrap(i - 1) |
|---|
| 285 | if (yb < min_lsp) then ! The minimum accepted is overtaken |
|---|
| 286 | nyears_maxlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 287 | found_lsp = .true. |
|---|
| 288 | else if (max_lsp < yb) then ! The maximum accepted is overtaken |
|---|
| 289 | nyears_maxlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year |
|---|
| 290 | found_lsp = .true. |
|---|
| 291 | endif |
|---|
| 292 | endif |
|---|
| 293 | |
|---|
| 294 | if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found |
|---|
| 295 | enddo |
|---|
| 296 | |
|---|
| 297 | deallocate(lsplask_unwrap) |
|---|
| 298 | |
|---|
| 299 | if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 300 | if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 301 | if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 302 | |
|---|
| 303 | write(*,*) 'Max. number of years for the obl. parameter =', nyears_maxobl |
|---|
| 304 | write(*,*) 'Max. number of years for the ecc. parameter =', nyears_maxecc |
|---|
| 305 | write(*,*) 'Max. number of years for the Lsp parameter =', nyears_maxlsp |
|---|
| 306 | |
|---|
| 307 | nyears_maxorb = min(nyears_maxobl,nyears_maxecc,nyears_maxlsp) |
|---|
| 308 | if (nyears_maxorb < 1.) nyears_maxorb = 1. ! nyears_maxorb should be at least equal to 1 |
|---|
| 309 | write(*,*) 'So the max. number of years for the orbital criteria =', nyears_maxorb |
|---|
| 310 | |
|---|
| 311 | END SUBROUTINE orbit_param_criterion |
|---|
| 312 | !======================================================================= |
|---|
| 313 | |
|---|
| 314 | !======================================================================= |
|---|
| 315 | SUBROUTINE set_new_orbitpm(i_myear,i_myear_leg) |
|---|
| 316 | !----------------------------------------------------------------------- |
|---|
| 317 | ! NAME |
|---|
| 318 | ! set_new_orbitpm |
|---|
| 319 | ! |
|---|
| 320 | ! DESCRIPTION |
|---|
| 321 | ! Recompute orbit parameters from Laskar values at a new year and |
|---|
| 322 | ! update planetary constants; handles Lsp modulo crossings. |
|---|
| 323 | ! |
|---|
| 324 | ! AUTHORS & DATE |
|---|
| 325 | ! R. Vandemeulebrouck |
|---|
| 326 | ! JB Clement, 2023-2025 |
|---|
| 327 | ! |
|---|
| 328 | ! NOTES |
|---|
| 329 | ! |
|---|
| 330 | !----------------------------------------------------------------------- |
|---|
| 331 | |
|---|
| 332 | use evolution, only: year_bp_ini, var_obl, var_ecc, var_lsp |
|---|
| 333 | use phys_constants, only: pi |
|---|
| 334 | use planete_h, only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day |
|---|
| 335 | use call_dayperi_mod, only: call_dayperi |
|---|
| 336 | |
|---|
| 337 | ! DECLARATION |
|---|
| 338 | ! ----------- |
|---|
| 339 | implicit none |
|---|
| 340 | |
|---|
| 341 | ! ARGUMENTS |
|---|
| 342 | ! --------- |
|---|
| 343 | real, intent(in) :: i_myear ! Number of simulated Martian years |
|---|
| 344 | real, intent(in) :: i_myear_leg ! Number of iterations of the PEM |
|---|
| 345 | |
|---|
| 346 | ! LOCAL VARIABLES |
|---|
| 347 | ! --------------- |
|---|
| 348 | real :: Year ! Year of the simulation |
|---|
| 349 | real :: Lsp ! Ls perihelion |
|---|
| 350 | real :: halfaxe ! Million of km |
|---|
| 351 | integer :: i |
|---|
| 352 | real :: unitastr |
|---|
| 353 | real :: xa, xb, ya, yb ! Variables for interpolation |
|---|
| 354 | logical :: found_year ! Flag variable |
|---|
| 355 | |
|---|
| 356 | ! CODE |
|---|
| 357 | ! ---- |
|---|
| 358 | ! ********************************************************************** |
|---|
| 359 | ! 0. Initializations |
|---|
| 360 | ! ********************************************************************** |
|---|
| 361 | Year = year_bp_ini + i_myear ! We compute the current year |
|---|
| 362 | Lsp = lsperi*180./pi ! We convert in degree |
|---|
| 363 | |
|---|
| 364 | write(*,*) 'set_new_orbitpm, Old year =', Year - i_myear_leg |
|---|
| 365 | write(*,*) 'set_new_orbitpm, Old obl. =', obliquit |
|---|
| 366 | write(*,*) 'set_new_orbitpm, Old ecc. =', e_elips |
|---|
| 367 | write(*,*) 'set_new_orbitpm, Old Lsp =', Lsp |
|---|
| 368 | write(*,*) 'set_new_orbitpm, New year =', Year |
|---|
| 369 | |
|---|
| 370 | found_year = .false. |
|---|
| 371 | if (Year < 0.) then |
|---|
| 372 | do i = iyear_lask,2,-1 |
|---|
| 373 | xa = year_lask(i) |
|---|
| 374 | xb = year_lask(i - 1) |
|---|
| 375 | if (xa <= Year .and. Year < xb) then |
|---|
| 376 | ! Obliquity |
|---|
| 377 | if (var_obl) then |
|---|
| 378 | ya = obl_lask(i) |
|---|
| 379 | yb = obl_lask(i - 1) |
|---|
| 380 | obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya |
|---|
| 381 | endif |
|---|
| 382 | |
|---|
| 383 | ! Eccentricity |
|---|
| 384 | if (var_ecc) then |
|---|
| 385 | ya = ecc_lask(i) |
|---|
| 386 | yb = ecc_lask(i - 1) |
|---|
| 387 | e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya |
|---|
| 388 | endif |
|---|
| 389 | |
|---|
| 390 | ! Lsp |
|---|
| 391 | if (var_lsp) then |
|---|
| 392 | ya = lsp_lask(i) |
|---|
| 393 | yb = lsp_lask(i - 1) |
|---|
| 394 | if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval |
|---|
| 395 | if (ya < yb) then ! Lsp should be decreasing |
|---|
| 396 | ya = ya + 360. |
|---|
| 397 | else ! Lsp should be increasing |
|---|
| 398 | yb = yb + 360. |
|---|
| 399 | endif |
|---|
| 400 | endif |
|---|
| 401 | Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.) |
|---|
| 402 | endif |
|---|
| 403 | found_year = .true. |
|---|
| 404 | exit ! The loop is left as soon as the right interval is found |
|---|
| 405 | endif |
|---|
| 406 | enddo |
|---|
| 407 | else if (Year >= 0 .and. Year < 100.) then ! For present orbital characteristics |
|---|
| 408 | if (var_obl) obliquit = 25.19 |
|---|
| 409 | if (var_ecc) e_elips = 0.0934 |
|---|
| 410 | if (var_lsp) Lsp = 251. |
|---|
| 411 | found_year = .true. |
|---|
| 412 | endif |
|---|
| 413 | |
|---|
| 414 | if (.not. found_year) error stop 'The new year could not be found in the file "obl_ecc_lsp.asc".' |
|---|
| 415 | |
|---|
| 416 | write(*,*) 'set_new_orbitpm, New obl. =', obliquit |
|---|
| 417 | write(*,*) 'set_new_orbitpm, New ecc. =', e_elips |
|---|
| 418 | write(*,*) 'set_new_orbitpm, New Lsp =', Lsp |
|---|
| 419 | |
|---|
| 420 | halfaxe = 227.94 |
|---|
| 421 | Lsp = Lsp*pi/180. |
|---|
| 422 | periheli = halfaxe*(1 - e_elips) |
|---|
| 423 | aphelie = halfaxe*(1 + e_elips) |
|---|
| 424 | unitastr = 149.597927 |
|---|
| 425 | p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr |
|---|
| 426 | |
|---|
| 427 | call call_dayperi(Lsp,e_elips,peri_day,year_day) |
|---|
| 428 | |
|---|
| 429 | call end_orbit() |
|---|
| 430 | |
|---|
| 431 | END SUBROUTINE set_new_orbitpm |
|---|
| 432 | !======================================================================= |
|---|
| 433 | |
|---|
| 434 | END MODULE orbit |
|---|
| 435 | |
|---|
