| 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 | ! DEPENDENCIES |
|---|
| 16 | ! ------------ |
|---|
| 17 | use numerics, only: dp, di, k4, largest_nb, minieps |
|---|
| 18 | |
|---|
| 19 | ! DECLARATION |
|---|
| 20 | ! ----------- |
|---|
| 21 | implicit none |
|---|
| 22 | |
|---|
| 23 | ! PARAMETERS |
|---|
| 24 | ! ---------- |
|---|
| 25 | ! File-related: |
|---|
| 26 | character(15), parameter :: orbitfile_name = 'obl_ecc_lsp.asc' |
|---|
| 27 | integer(di), protected :: iyear_lask! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year |
|---|
| 28 | real(dp), dimension(:), allocatable, protected :: year_lask ! Current_year before present from Laskar et al. |
|---|
| 29 | real(dp), dimension(:), allocatable, protected :: obl_lask ! Obliquity in Laskar's data [deg] |
|---|
| 30 | real(dp), dimension(:), allocatable, protected :: ecc_lask ! Eccentricity in Laskar's data |
|---|
| 31 | real(dp), dimension(:), allocatable, protected :: lsp_lask ! Ls perihelie in Laskar's data [deg] |
|---|
| 32 | ! Planet-related: |
|---|
| 33 | real(dp), protected :: obliquity ! Planet obliquity [deg] |
|---|
| 34 | real(dp), protected :: eccentricity ! Orbit eccentricity |
|---|
| 35 | real(dp), protected :: ls_peri ! Solar longitude of the perihelion [deg] |
|---|
| 36 | real(dp), protected :: perihelion ! Perihelion [Mkm] |
|---|
| 37 | real(dp), protected :: semimaj_axis ! Semi-major axis [Mkm] |
|---|
| 38 | real(dp), protected :: aphelion ! Aphelion [Mkm] |
|---|
| 39 | real(dp), protected :: date_peri ! Date of perihelion [sol] |
|---|
| 40 | real(dp), protected :: sol_len_s ! Length of a sol [s] |
|---|
| 41 | real(dp), protected :: yr_len_sol ! Length of a year [sol] |
|---|
| 42 | ! Evolution-related: |
|---|
| 43 | logical(k4), protected :: evol_orbit ! Flag to follow the orbital parameters |
|---|
| 44 | logical(k4), protected :: evol_obl ! Flag the PEM to follow obliquity |
|---|
| 45 | logical(k4), protected :: evol_ecc ! Flag the PEM to follow eccenticity |
|---|
| 46 | logical(k4), protected :: evol_lsp ! Flag the PEM to follow Ls perihelie |
|---|
| 47 | real(dp), protected :: max_change_obl ! Maximum admissible change for obliquity |
|---|
| 48 | real(dp), protected :: max_change_ecc ! Maximum admissible change for eccentricity |
|---|
| 49 | real(dp), protected :: max_change_lsp ! Maximum admissible change for Lsp |
|---|
| 50 | |
|---|
| 51 | contains |
|---|
| 52 | !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|---|
| 53 | |
|---|
| 54 | !======================================================================= |
|---|
| 55 | SUBROUTINE set_orbit_config(evol_orbit_in,evol_obl_in,evol_ecc_in,evol_lsp_in,max_change_obl_in,max_change_ecc_in,max_change_lsp_in) |
|---|
| 56 | !----------------------------------------------------------------------- |
|---|
| 57 | ! NAME |
|---|
| 58 | ! set_orbit_config |
|---|
| 59 | ! |
|---|
| 60 | ! DESCRIPTION |
|---|
| 61 | ! Setter for 'orbit' configuration parameters. |
|---|
| 62 | ! |
|---|
| 63 | ! AUTHORS & DATE |
|---|
| 64 | ! JB Clement, 02/2026 |
|---|
| 65 | ! |
|---|
| 66 | ! NOTES |
|---|
| 67 | ! |
|---|
| 68 | !----------------------------------------------------------------------- |
|---|
| 69 | |
|---|
| 70 | ! DEPENDENCIES |
|---|
| 71 | ! ------------ |
|---|
| 72 | use utility, only: real2str, bool2str |
|---|
| 73 | use display, only: print_msg |
|---|
| 74 | |
|---|
| 75 | ! DECLARATION |
|---|
| 76 | ! ----------- |
|---|
| 77 | implicit none |
|---|
| 78 | |
|---|
| 79 | ! ARGUMENTS |
|---|
| 80 | ! --------- |
|---|
| 81 | logical(k4), intent(in) :: evol_orbit_in, evol_obl_in, evol_ecc_in, evol_lsp_in |
|---|
| 82 | real(dp), intent(in) :: max_change_obl_in, max_change_ecc_in, max_change_lsp_in |
|---|
| 83 | |
|---|
| 84 | ! CODE |
|---|
| 85 | ! ---- |
|---|
| 86 | evol_orbit = evol_orbit_in |
|---|
| 87 | evol_obl = evol_obl_in |
|---|
| 88 | evol_ecc = evol_ecc_in |
|---|
| 89 | evol_lsp = evol_lsp_in |
|---|
| 90 | max_change_obl = max_change_obl_in |
|---|
| 91 | max_change_ecc = max_change_ecc_in |
|---|
| 92 | max_change_lsp = max_change_lsp_in |
|---|
| 93 | call print_msg('evol_orbit = '//bool2str(evol_orbit)) |
|---|
| 94 | call print_msg('evol_obl = '//bool2str(evol_obl)) |
|---|
| 95 | call print_msg('evol_ecc = '//bool2str(evol_ecc)) |
|---|
| 96 | call print_msg('evol_lsp = '//bool2str(evol_lsp)) |
|---|
| 97 | call print_msg('max_change_obl = '//real2str(max_change_obl)) |
|---|
| 98 | call print_msg('max_change_ecc = '//real2str(max_change_ecc)) |
|---|
| 99 | call print_msg('max_change_lsp = '//real2str(max_change_lsp)) |
|---|
| 100 | |
|---|
| 101 | END SUBROUTINE set_orbit_config |
|---|
| 102 | !======================================================================= |
|---|
| 103 | |
|---|
| 104 | !======================================================================= |
|---|
| 105 | SUBROUTINE init_orbit(obliquity_in,perihelion_in,aphelion_in,date_peri_in,yr_len_sol_in,sol_len_s_in) |
|---|
| 106 | !----------------------------------------------------------------------- |
|---|
| 107 | ! NAME |
|---|
| 108 | ! init_orbit |
|---|
| 109 | ! |
|---|
| 110 | ! DESCRIPTION |
|---|
| 111 | ! Initialize orbital characteristics of the planet. |
|---|
| 112 | ! |
|---|
| 113 | ! AUTHORS & DATE |
|---|
| 114 | ! JB Clement, 12/2025 |
|---|
| 115 | ! |
|---|
| 116 | ! NOTES |
|---|
| 117 | ! Taken from 'iniorbit' in the Mars PCM. |
|---|
| 118 | !----------------------------------------------------------------------- |
|---|
| 119 | |
|---|
| 120 | ! DEPENDENCIES |
|---|
| 121 | ! ------------ |
|---|
| 122 | use maths, only: pi |
|---|
| 123 | use display, only: print_msg |
|---|
| 124 | use utility, only: real2str |
|---|
| 125 | |
|---|
| 126 | ! DECLARATION |
|---|
| 127 | ! ----------- |
|---|
| 128 | implicit none |
|---|
| 129 | |
|---|
| 130 | ! ARGUMENTS |
|---|
| 131 | ! --------- |
|---|
| 132 | real(dp), intent(in) :: obliquity_in, perihelion_in, aphelion_in, date_peri_in, yr_len_sol_in, sol_len_s_in |
|---|
| 133 | |
|---|
| 134 | ! LOCAL VARIABLES |
|---|
| 135 | ! --------------- |
|---|
| 136 | real(dp) :: zxref, zanom, zz, zx0, zdx |
|---|
| 137 | real(dp) :: lsp ! Solar longitude of the perihelion [rad] |
|---|
| 138 | integer(di) :: iter |
|---|
| 139 | |
|---|
| 140 | ! CODE |
|---|
| 141 | ! ---- |
|---|
| 142 | ! Setters |
|---|
| 143 | obliquity = obliquity_in |
|---|
| 144 | perihelion = perihelion_in |
|---|
| 145 | aphelion = aphelion_in |
|---|
| 146 | date_peri = date_peri_in |
|---|
| 147 | sol_len_s = sol_len_s_in |
|---|
| 148 | yr_len_sol = yr_len_sol_in |
|---|
| 149 | |
|---|
| 150 | ! Compute eccentricity and semi-major axis |
|---|
| 151 | eccentricity = (aphelion - perihelion)/(aphelion + perihelion) |
|---|
| 152 | semimaj_axis = 0.5_dp*(aphelion + perihelion) |
|---|
| 153 | |
|---|
| 154 | ! Compute mean anomaly zanom |
|---|
| 155 | zz = (yr_len_sol - date_peri)/yr_len_sol |
|---|
| 156 | zanom = 2._dp*pi*(zz - nint(zz)) |
|---|
| 157 | zxref = abs(zanom) |
|---|
| 158 | |
|---|
| 159 | ! Solve equation zx0 - e*sin(zx0) = zxref for eccentric anomaly zx0 using Newton method |
|---|
| 160 | zx0 = zxref + eccentricity*sin(zxref) |
|---|
| 161 | do iter = 1,100 |
|---|
| 162 | zdx = -(zx0 - eccentricity*sin(zx0) - zxref)/(1._dp - eccentricity*cos(zx0)) |
|---|
| 163 | if (abs(zdx) <= 1.e-12_dp) exit |
|---|
| 164 | zx0 = zx0 + zdx |
|---|
| 165 | end do |
|---|
| 166 | |
|---|
| 167 | zx0 = zx0 + zdx |
|---|
| 168 | if (zanom < 0.) zx0 = -zx0 |
|---|
| 169 | |
|---|
| 170 | ! Compute solar longitude of the perihelion |
|---|
| 171 | lsp = -2._dp*atan(sqrt((1._dp + eccentricity)/(1._dp - eccentricity))*tan(zx0/2._dp)) |
|---|
| 172 | lsp = modulo(lsp,2._dp*pi) |
|---|
| 173 | ls_peri = lsp*180._dp/pi ! Conversion in degree |
|---|
| 174 | |
|---|
| 175 | call print_msg('Obliquity [deg] = '//real2str(obliquity)) |
|---|
| 176 | call print_msg('Eccentricity [-] = '//real2str(eccentricity)) |
|---|
| 177 | call print_msg('Ls of perihelion [deg] = '//real2str(ls_peri)) |
|---|
| 178 | call print_msg('Perihelion [Mkm] = '//real2str(perihelion)) |
|---|
| 179 | call print_msg('Aphelion [Mkm] = '//real2str(aphelion)) |
|---|
| 180 | call print_msg('Semi-major axis [Mkm] = '//real2str(semimaj_axis)) |
|---|
| 181 | call print_msg('Date of perihelion [sol] = '//real2str(date_peri)) |
|---|
| 182 | call print_msg('Length of a sol [s] = '//real2str(sol_len_s)) |
|---|
| 183 | call print_msg('Length of a year [sol] = '//real2str(yr_len_sol)) |
|---|
| 184 | |
|---|
| 185 | END SUBROUTINE init_orbit |
|---|
| 186 | !======================================================================= |
|---|
| 187 | |
|---|
| 188 | !======================================================================= |
|---|
| 189 | SUBROUTINE lsp2datep(lsp,datep) |
|---|
| 190 | !----------------------------------------------------------------------- |
|---|
| 191 | ! NAME |
|---|
| 192 | ! lsp2datep |
|---|
| 193 | ! |
|---|
| 194 | ! DESCRIPTION |
|---|
| 195 | ! Initialize orbital data of the planet. |
|---|
| 196 | ! |
|---|
| 197 | ! AUTHORS & DATE |
|---|
| 198 | ! JB Clement, 12/2025 |
|---|
| 199 | ! |
|---|
| 200 | ! NOTES |
|---|
| 201 | ! Taken from 'lsp2solp' in the Mars PCM. |
|---|
| 202 | !----------------------------------------------------------------------- |
|---|
| 203 | |
|---|
| 204 | ! DEPENDENCIES |
|---|
| 205 | ! ------------ |
|---|
| 206 | use maths, only: pi |
|---|
| 207 | |
|---|
| 208 | ! DECLARATION |
|---|
| 209 | ! ----------- |
|---|
| 210 | implicit none |
|---|
| 211 | |
|---|
| 212 | ! ARGUMENTS |
|---|
| 213 | ! --------- |
|---|
| 214 | real(dp), intent(in) :: lsp ! Ls at perihelion [deg] |
|---|
| 215 | real(dp), intent(out) :: datep ! Date at perihelion [sol] |
|---|
| 216 | |
|---|
| 217 | ! LOCAL VARIABLES |
|---|
| 218 | ! --------------- |
|---|
| 219 | real(dp) :: zx0 ! Eccentric anomaly at Ls = 0 |
|---|
| 220 | |
|---|
| 221 | ! CODE |
|---|
| 222 | ! ---- |
|---|
| 223 | ! Compute orbit geometrical parameters |
|---|
| 224 | zx0 = -2._dp*atan(tan(0.5_dp*lsp*pi/180._dp)*sqrt((1._dp - eccentricity)/(1._dp + eccentricity))) |
|---|
| 225 | if (zx0 <= 0._dp) zx0 = zx0 + 2._dp*pi |
|---|
| 226 | |
|---|
| 227 | ! Compute sol at perihelion |
|---|
| 228 | datep = yr_len_sol*(1._dp - (zx0 - eccentricity*sin(zx0))/(2._dp*pi)) |
|---|
| 229 | |
|---|
| 230 | END SUBROUTINE lsp2datep |
|---|
| 231 | !======================================================================= |
|---|
| 232 | |
|---|
| 233 | !======================================================================= |
|---|
| 234 | SUBROUTINE ini_orbit() |
|---|
| 235 | !----------------------------------------------------------------------- |
|---|
| 236 | ! NAME |
|---|
| 237 | ! ini_orbit |
|---|
| 238 | ! |
|---|
| 239 | ! DESCRIPTION |
|---|
| 240 | ! Initialize the parameters of module 'orbit'. |
|---|
| 241 | ! |
|---|
| 242 | ! AUTHORS & DATE |
|---|
| 243 | ! R. Vandemeulebrouck |
|---|
| 244 | ! JB Clement, 2023-2025 |
|---|
| 245 | ! |
|---|
| 246 | ! NOTES |
|---|
| 247 | ! |
|---|
| 248 | !----------------------------------------------------------------------- |
|---|
| 249 | |
|---|
| 250 | ! DEPENDENCIES |
|---|
| 251 | ! ------------ |
|---|
| 252 | use stoppage, only: stop_clean |
|---|
| 253 | |
|---|
| 254 | ! DECLARATION |
|---|
| 255 | ! ----------- |
|---|
| 256 | implicit none |
|---|
| 257 | |
|---|
| 258 | ! LOCAL VARIABLES |
|---|
| 259 | ! --------------- |
|---|
| 260 | integer(di) :: n ! Number of lines in the file |
|---|
| 261 | integer(di) :: ierr, funit |
|---|
| 262 | logical(k4) :: here |
|---|
| 263 | |
|---|
| 264 | ! CODE |
|---|
| 265 | ! ---- |
|---|
| 266 | inquire(file = orbitfile_name,exist = here) |
|---|
| 267 | if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//orbitfile_name//'"!',1) |
|---|
| 268 | |
|---|
| 269 | ! Get the number of lines for Laskar's orbital file |
|---|
| 270 | open(newunit = funit,file = orbitfile_name,status = 'old',action = 'read',iostat = ierr) |
|---|
| 271 | if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//orbitfile_name//'"!',ierr) |
|---|
| 272 | n = 0 |
|---|
| 273 | do |
|---|
| 274 | read(funit,*,iostat = ierr) |
|---|
| 275 | if (ierr /= 0) exit |
|---|
| 276 | n = n + 1 |
|---|
| 277 | end do |
|---|
| 278 | close(funit) |
|---|
| 279 | |
|---|
| 280 | ! Allocation of module arrays |
|---|
| 281 | if (.not. allocated(year_lask)) allocate(year_lask(n)) |
|---|
| 282 | if (.not. allocated(obl_lask)) allocate(obl_lask(n)) |
|---|
| 283 | if (.not. allocated(ecc_lask)) allocate(ecc_lask(n)) |
|---|
| 284 | if (.not. allocated(lsp_lask)) allocate(lsp_lask(n)) |
|---|
| 285 | year_lask(:) = 0._dp |
|---|
| 286 | obl_lask(:) = 0._dp |
|---|
| 287 | ecc_lask(:) = 0._dp |
|---|
| 288 | lsp_lask(:) = 0._dp |
|---|
| 289 | |
|---|
| 290 | END SUBROUTINE ini_orbit |
|---|
| 291 | !======================================================================= |
|---|
| 292 | |
|---|
| 293 | !======================================================================= |
|---|
| 294 | SUBROUTINE end_orbit |
|---|
| 295 | !----------------------------------------------------------------------- |
|---|
| 296 | ! NAME |
|---|
| 297 | ! end_orbit |
|---|
| 298 | ! |
|---|
| 299 | ! DESCRIPTION |
|---|
| 300 | ! Deallocate orbital data arrays. |
|---|
| 301 | ! |
|---|
| 302 | ! AUTHORS & DATE |
|---|
| 303 | ! R. Vandemeulebrouck |
|---|
| 304 | ! JB Clement, 2023-2025 |
|---|
| 305 | ! |
|---|
| 306 | ! NOTES |
|---|
| 307 | ! |
|---|
| 308 | !----------------------------------------------------------------------- |
|---|
| 309 | |
|---|
| 310 | ! DECLARATION |
|---|
| 311 | ! ----------- |
|---|
| 312 | implicit none |
|---|
| 313 | |
|---|
| 314 | ! CODE |
|---|
| 315 | ! ---- |
|---|
| 316 | if (allocated(year_lask)) deallocate(year_lask) |
|---|
| 317 | if (allocated(obl_lask)) deallocate(obl_lask) |
|---|
| 318 | if (allocated(ecc_lask)) deallocate(ecc_lask) |
|---|
| 319 | if (allocated(lsp_lask)) deallocate(lsp_lask) |
|---|
| 320 | |
|---|
| 321 | END SUBROUTINE end_orbit |
|---|
| 322 | !======================================================================= |
|---|
| 323 | |
|---|
| 324 | !======================================================================= |
|---|
| 325 | SUBROUTINE read_orbitpm(n_yr_sim) |
|---|
| 326 | !----------------------------------------------------------------------- |
|---|
| 327 | ! NAME |
|---|
| 328 | ! read_orbitpm |
|---|
| 329 | ! |
|---|
| 330 | ! DESCRIPTION |
|---|
| 331 | ! Read Laskar's orbital file and locate index for |
|---|
| 332 | ! closest lower year relative to current simulation year. |
|---|
| 333 | ! |
|---|
| 334 | ! AUTHORS & DATE |
|---|
| 335 | ! JB Clement, 2023-2025 |
|---|
| 336 | ! |
|---|
| 337 | ! NOTES |
|---|
| 338 | ! |
|---|
| 339 | !----------------------------------------------------------------------- |
|---|
| 340 | |
|---|
| 341 | ! DEPENDENCIES |
|---|
| 342 | ! ------------ |
|---|
| 343 | use evolution, only: pem_ini_date, r_plnt2earth_yr |
|---|
| 344 | use stoppage, only: stop_clean |
|---|
| 345 | use display, only: print_msg |
|---|
| 346 | use utility, only: real2str, int2str |
|---|
| 347 | |
|---|
| 348 | ! DECLARATION |
|---|
| 349 | ! ----------- |
|---|
| 350 | implicit none |
|---|
| 351 | |
|---|
| 352 | ! ARGUMENTS |
|---|
| 353 | ! --------- |
|---|
| 354 | real(dp), intent(in) :: n_yr_sim |
|---|
| 355 | |
|---|
| 356 | ! LOCAL VARIABLES |
|---|
| 357 | ! --------------- |
|---|
| 358 | integer(di) :: ierr, funit, i, n |
|---|
| 359 | logical(k4) :: here |
|---|
| 360 | real(dp) :: current_year |
|---|
| 361 | |
|---|
| 362 | ! CODE |
|---|
| 363 | ! ---- |
|---|
| 364 | call print_msg('> Reading "'//orbitfile_name//'"') |
|---|
| 365 | ! Compute the current year |
|---|
| 366 | current_year = pem_ini_date + n_yr_sim |
|---|
| 367 | call print_msg('Current year = '//real2str(current_year)) |
|---|
| 368 | |
|---|
| 369 | inquire(file = orbitfile_name,exist = here) |
|---|
| 370 | if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//orbitfile_name//'"!',1) |
|---|
| 371 | |
|---|
| 372 | ! Read the Laskars's orbital file |
|---|
| 373 | open(newunit = funit,file = orbitfile_name,status = 'old',action = 'read',iostat = ierr) |
|---|
| 374 | if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//orbitfile_name//'"!',ierr) |
|---|
| 375 | n = size(year_lask,1) ! Number of lines in the file |
|---|
| 376 | iyear_lask = 0 ! Line number for closest lower year relative to current simulation year |
|---|
| 377 | do i = 1,n |
|---|
| 378 | read(funit,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i) |
|---|
| 379 | year_lask(i) = year_lask(i)*1000._dp/r_plnt2earth_yr ! Conversion from Laskar's kilo Earth years to planetary years |
|---|
| 380 | if (year_lask(i) > current_year) iyear_lask = i + 1 |
|---|
| 381 | end do |
|---|
| 382 | close(funit) |
|---|
| 383 | if (iyear_lask == 0 .or. iyear_lask == n + 1) then |
|---|
| 384 | call stop_clean(__FILE__,__LINE__,'the current year could not be found in the file "'//orbitfile_name//'".',1) |
|---|
| 385 | else |
|---|
| 386 | call print_msg('The current year in "'//orbitfile_name//'" is at line '//int2str(iyear_lask)//'.') |
|---|
| 387 | end if |
|---|
| 388 | |
|---|
| 389 | END SUBROUTINE read_orbitpm |
|---|
| 390 | !======================================================================= |
|---|
| 391 | |
|---|
| 392 | !======================================================================= |
|---|
| 393 | SUBROUTINE compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb) |
|---|
| 394 | !----------------------------------------------------------------------- |
|---|
| 395 | ! NAME |
|---|
| 396 | ! compute_maxyr_orbit |
|---|
| 397 | ! |
|---|
| 398 | ! DESCRIPTION |
|---|
| 399 | ! Determine maximum number of years before orbital params change |
|---|
| 400 | ! beyond admissible thresholds. |
|---|
| 401 | ! |
|---|
| 402 | ! AUTHORS & DATE |
|---|
| 403 | ! R. Vandemeulebrouck |
|---|
| 404 | ! JB Clement, 2023-2025 |
|---|
| 405 | ! |
|---|
| 406 | ! NOTES |
|---|
| 407 | ! Handles Lsp modulo crossings. |
|---|
| 408 | !----------------------------------------------------------------------- |
|---|
| 409 | |
|---|
| 410 | ! DEPENDENCIES |
|---|
| 411 | ! ------------ |
|---|
| 412 | use evolution, only: pem_ini_date |
|---|
| 413 | use display, only: print_msg |
|---|
| 414 | use utility, only: real2str |
|---|
| 415 | |
|---|
| 416 | ! DECLARATION |
|---|
| 417 | ! ----------- |
|---|
| 418 | implicit none |
|---|
| 419 | |
|---|
| 420 | ! ARGUMENTS |
|---|
| 421 | ! --------- |
|---|
| 422 | real(dp), intent(in) :: n_yr_sim ! Current year of the simulation |
|---|
| 423 | real(dp), intent(out) :: nmax_yr_runorb ! Maximum number of years of the PEM before orbital parameters change too much |
|---|
| 424 | |
|---|
| 425 | ! LOCAL VARIABLES |
|---|
| 426 | ! --------------- |
|---|
| 427 | real(dp) :: current_year |
|---|
| 428 | real(dp) :: max_obl, max_ecc, max_lsp |
|---|
| 429 | real(dp) :: min_obl, min_ecc, min_lsp |
|---|
| 430 | real(dp) :: nmax_yr_runobl, nmax_yr_runecc, nmax_yr_runlsp ! Maximum year iteration before reaching an inadmissible value of orbit param |
|---|
| 431 | real(dp) :: xa, xb, ya, yb |
|---|
| 432 | logical(k4) :: found_obl, found_ecc, found_lsp |
|---|
| 433 | real(dp) :: lsp_corr ! Lsp corrected if the "modulo is crossed" |
|---|
| 434 | real(dp) :: delta |
|---|
| 435 | real(dp), dimension(:), allocatable :: lsplask_unwrap |
|---|
| 436 | integer(di) :: i |
|---|
| 437 | |
|---|
| 438 | ! CODE |
|---|
| 439 | ! ---- |
|---|
| 440 | if (.not. evol_orbit) then |
|---|
| 441 | nmax_yr_runorb = largest_nb ! Default huge value |
|---|
| 442 | return |
|---|
| 443 | end if |
|---|
| 444 | call print_msg('> Computing the accepted maximum number of years due to orbital parameters') |
|---|
| 445 | current_year = pem_ini_date + n_yr_sim ! We compute the current year |
|---|
| 446 | |
|---|
| 447 | ! Unwrap the Lsp changes in Laskar's file |
|---|
| 448 | allocate(lsplask_unwrap(iyear_lask)) |
|---|
| 449 | lsplask_unwrap(iyear_lask) = lsp_lask(iyear_lask) |
|---|
| 450 | do i = iyear_lask,2,-1 |
|---|
| 451 | delta = lsp_lask(i - 1) - lsp_lask(i) |
|---|
| 452 | if (delta > 180._dp) then |
|---|
| 453 | lsp_corr = lsp_lask(i - 1) - 360._dp |
|---|
| 454 | else if (delta < -180._dp) then |
|---|
| 455 | lsp_corr = lsp_lask(i - 1) + 360._dp |
|---|
| 456 | else |
|---|
| 457 | lsp_corr = lsp_lask(i - 1) |
|---|
| 458 | end if |
|---|
| 459 | lsplask_unwrap(i - 1) = lsplask_unwrap(i) + lsp_corr - lsp_lask(i) |
|---|
| 460 | end do |
|---|
| 461 | |
|---|
| 462 | ! Correct the current Lsp according to the unwrapping |
|---|
| 463 | delta = ls_peri - lsplask_unwrap(iyear_lask) |
|---|
| 464 | if (delta > 180._dp) then |
|---|
| 465 | ls_peri = ls_peri - 360._dp |
|---|
| 466 | else if (delta < -180._dp) then |
|---|
| 467 | ls_peri = ls_peri + 360._dp |
|---|
| 468 | end if |
|---|
| 469 | |
|---|
| 470 | ! Maximum change accepted for orbital parameters |
|---|
| 471 | max_obl = obliquity + max_change_obl |
|---|
| 472 | min_obl = obliquity - max_change_obl |
|---|
| 473 | call print_msg('Obl. (current|accepted min|accepted max): '//real2str(obliquity)//' | '//real2str(min_obl)//real2str(max_obl)) |
|---|
| 474 | |
|---|
| 475 | max_ecc = min(eccentricity + max_change_ecc,1. - minieps) ! ecc < 1. |
|---|
| 476 | min_ecc = max(eccentricity - max_change_ecc,0._dp) ! ecc >= 0. |
|---|
| 477 | call print_msg('Ecc. (current|accepted min|accepted max): '//real2str(eccentricity)//' | '//real2str(min_ecc)//' | '//real2str(max_ecc)) |
|---|
| 478 | |
|---|
| 479 | max_lsp = ls_peri + max_change_lsp |
|---|
| 480 | min_lsp = ls_peri - max_change_lsp |
|---|
| 481 | call print_msg('Lsp (current|accepted min|accepted max): '//real2str(ls_peri)//' | '//real2str(min_lsp)//' | '//real2str(max_lsp)) |
|---|
| 482 | |
|---|
| 483 | ! Initialization |
|---|
| 484 | nmax_yr_runobl = largest_nb ! Default huge value |
|---|
| 485 | nmax_yr_runecc = largest_nb ! Default huge value |
|---|
| 486 | nmax_yr_runlsp =largest_nb ! Default huge value |
|---|
| 487 | found_obl = .false. |
|---|
| 488 | found_ecc = .false. |
|---|
| 489 | found_lsp = .false. |
|---|
| 490 | if (.not. evol_obl) found_obl = .true. |
|---|
| 491 | if (.not. evol_ecc) found_ecc = .true. |
|---|
| 492 | if (.not. evol_lsp) found_lsp = .true. |
|---|
| 493 | |
|---|
| 494 | ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters |
|---|
| 495 | do i = iyear_lask,2,-1 |
|---|
| 496 | xa = year_lask(i) |
|---|
| 497 | xb = year_lask(i - 1) |
|---|
| 498 | |
|---|
| 499 | ! Obliquity |
|---|
| 500 | if (evol_obl .and. .not. found_obl) then |
|---|
| 501 | ya = obl_lask(i) |
|---|
| 502 | yb = obl_lask(i - 1) |
|---|
| 503 | if (yb < min_obl) then ! The minimum accepted is overtaken |
|---|
| 504 | nmax_yr_runobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - current_year |
|---|
| 505 | found_obl = .true. |
|---|
| 506 | else if (max_obl < yb) then ! The maximum accepted is overtaken |
|---|
| 507 | nmax_yr_runobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - current_year |
|---|
| 508 | found_obl = .true. |
|---|
| 509 | end if |
|---|
| 510 | end if |
|---|
| 511 | |
|---|
| 512 | ! Eccentricity |
|---|
| 513 | if (evol_ecc .and. .not. found_ecc) then |
|---|
| 514 | ya = ecc_lask(i) |
|---|
| 515 | yb = ecc_lask(i - 1) |
|---|
| 516 | if (yb < min_ecc) then ! The minimum accepted is overtaken |
|---|
| 517 | nmax_yr_runecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - current_year |
|---|
| 518 | found_ecc = .true. |
|---|
| 519 | else if (max_ecc < yb) then ! The maximum accepted is overtaken |
|---|
| 520 | nmax_yr_runecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - current_year |
|---|
| 521 | found_ecc = .true. |
|---|
| 522 | end if |
|---|
| 523 | end if |
|---|
| 524 | |
|---|
| 525 | ! Lsp |
|---|
| 526 | if (evol_lsp .and. .not. found_lsp) then |
|---|
| 527 | ya = lsplask_unwrap(i) |
|---|
| 528 | yb = lsplask_unwrap(i - 1) |
|---|
| 529 | if (yb < min_lsp) then ! The minimum accepted is overtaken |
|---|
| 530 | nmax_yr_runlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - current_year |
|---|
| 531 | found_lsp = .true. |
|---|
| 532 | else if (max_lsp < yb) then ! The maximum accepted is overtaken |
|---|
| 533 | nmax_yr_runlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - current_year |
|---|
| 534 | found_lsp = .true. |
|---|
| 535 | end if |
|---|
| 536 | end if |
|---|
| 537 | |
|---|
| 538 | if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found |
|---|
| 539 | end do |
|---|
| 540 | |
|---|
| 541 | deallocate(lsplask_unwrap) |
|---|
| 542 | |
|---|
| 543 | if (.not. found_obl) call print_msg('The maximum reachable year for obl. could not be found in "'//orbitfile_name//'".') |
|---|
| 544 | if (.not. found_ecc) call print_msg('The maximum reachable year for ecc. could not be found in "'//orbitfile_name//'".') |
|---|
| 545 | if (.not. found_lsp) call print_msg('The maximum reachable year for Lsp could not be found in "'//orbitfile_name//'".') |
|---|
| 546 | |
|---|
| 547 | call print_msg('Number of years accepted for obl. = '//real2str(nmax_yr_runobl)) |
|---|
| 548 | call print_msg('Number of years accepted for ecc. = '//real2str(nmax_yr_runecc)) |
|---|
| 549 | call print_msg('Number of years accepted for Lsp = '//real2str(nmax_yr_runlsp)) |
|---|
| 550 | |
|---|
| 551 | nmax_yr_runorb = min(nmax_yr_runobl,nmax_yr_runecc,nmax_yr_runlsp) |
|---|
| 552 | if (nmax_yr_runorb < 1._dp) nmax_yr_runorb = 1._dp ! nmax_yr_runorb should be at least equal to 1 |
|---|
| 553 | |
|---|
| 554 | call print_msg('So the max. number of years for the orbital criterion = '//real2str(nmax_yr_runorb)) |
|---|
| 555 | |
|---|
| 556 | END SUBROUTINE compute_maxyr_orbit |
|---|
| 557 | !======================================================================= |
|---|
| 558 | |
|---|
| 559 | !======================================================================= |
|---|
| 560 | SUBROUTINE update_orbit(n_yr_sim,n_yr_run) |
|---|
| 561 | !----------------------------------------------------------------------- |
|---|
| 562 | ! NAME |
|---|
| 563 | ! update_orbit |
|---|
| 564 | ! |
|---|
| 565 | ! DESCRIPTION |
|---|
| 566 | ! Update orbital parameters from Laskar values at a new year. |
|---|
| 567 | ! |
|---|
| 568 | ! AUTHORS & DATE |
|---|
| 569 | ! R. Vandemeulebrouck |
|---|
| 570 | ! JB Clement, 2023-2025 |
|---|
| 571 | ! |
|---|
| 572 | ! NOTES |
|---|
| 573 | ! |
|---|
| 574 | !----------------------------------------------------------------------- |
|---|
| 575 | |
|---|
| 576 | use evolution, only: pem_ini_date |
|---|
| 577 | use call_dayperi_mod, only: call_dayperi |
|---|
| 578 | use stoppage, only: stop_clean |
|---|
| 579 | use display, only: print_msg |
|---|
| 580 | use utility, only: real2str |
|---|
| 581 | |
|---|
| 582 | ! DECLARATION |
|---|
| 583 | ! ----------- |
|---|
| 584 | implicit none |
|---|
| 585 | |
|---|
| 586 | ! ARGUMENTS |
|---|
| 587 | ! --------- |
|---|
| 588 | real(dp), intent(in) :: n_yr_sim ! Number of simulated years |
|---|
| 589 | real(dp), intent(in) :: n_yr_run ! Number of iterations of the PEM |
|---|
| 590 | |
|---|
| 591 | ! LOCAL VARIABLES |
|---|
| 592 | ! --------------- |
|---|
| 593 | integer(di) :: i |
|---|
| 594 | real(dp) :: current_year, old_year, obl_old, ecc_old, lsp_old |
|---|
| 595 | real(dp) :: xa, xb, ya, yb ! Variables for interpolation |
|---|
| 596 | logical(k4) :: found_year |
|---|
| 597 | |
|---|
| 598 | ! CODE |
|---|
| 599 | ! ---- |
|---|
| 600 | call print_msg('> Updating the orbital parameters for the new year') |
|---|
| 601 | current_year = pem_ini_date + n_yr_sim ! We compute the current year |
|---|
| 602 | old_year = current_year - n_yr_run |
|---|
| 603 | obl_old = obliquity |
|---|
| 604 | ecc_old = eccentricity |
|---|
| 605 | lsp_old = ls_peri |
|---|
| 606 | |
|---|
| 607 | found_year = .false. |
|---|
| 608 | if (current_year < 0._dp) then |
|---|
| 609 | do i = iyear_lask,2,-1 |
|---|
| 610 | xa = year_lask(i) |
|---|
| 611 | xb = year_lask(i - 1) |
|---|
| 612 | if (xa <= current_year .and. current_year < xb) then |
|---|
| 613 | ! Obliquity |
|---|
| 614 | if (evol_obl) then |
|---|
| 615 | ya = obl_lask(i) |
|---|
| 616 | yb = obl_lask(i - 1) |
|---|
| 617 | obliquity = (current_year - xa)*(yb - ya)/(xb - xa) + ya |
|---|
| 618 | end if |
|---|
| 619 | |
|---|
| 620 | ! Eccentricity |
|---|
| 621 | if (evol_ecc) then |
|---|
| 622 | ya = ecc_lask(i) |
|---|
| 623 | yb = ecc_lask(i - 1) |
|---|
| 624 | eccentricity = (current_year - xa)*(yb - ya)/(xb - xa) + ya |
|---|
| 625 | end if |
|---|
| 626 | |
|---|
| 627 | ! Lsp |
|---|
| 628 | if (evol_lsp) then |
|---|
| 629 | ya = lsp_lask(i) |
|---|
| 630 | yb = lsp_lask(i - 1) |
|---|
| 631 | if (abs(yb - ya) > 300._dp) then ! If modulo is "crossed" through the interval |
|---|
| 632 | if (ya < yb) then ! Lsp should be decreasing |
|---|
| 633 | ya = ya + 360._dp |
|---|
| 634 | else ! Lsp should be increasing |
|---|
| 635 | yb = yb + 360._dp |
|---|
| 636 | end if |
|---|
| 637 | end if |
|---|
| 638 | ls_peri = modulo((current_year - xa)*(yb - ya)/(xb - xa) + ya,360._dp) |
|---|
| 639 | end if |
|---|
| 640 | found_year = .true. |
|---|
| 641 | exit ! The loop is left as soon as the right interval is found |
|---|
| 642 | end if |
|---|
| 643 | end do |
|---|
| 644 | else if (current_year >= 0._dp .and. current_year < 100._dp) then ! For present orbital characteristics |
|---|
| 645 | if (evol_obl) obliquity = obl_lask(1) |
|---|
| 646 | if (evol_ecc) eccentricity = ecc_lask(1) |
|---|
| 647 | if (evol_lsp) ls_peri = lsp_lask(1) |
|---|
| 648 | found_year = .true. |
|---|
| 649 | end if |
|---|
| 650 | |
|---|
| 651 | if (.not. found_year) call stop_clean(__FILE__,__LINE__,'the new year could not be found in the file "'//orbitfile_name//'".',1) |
|---|
| 652 | |
|---|
| 653 | call lsp2datep(ls_peri,date_peri) |
|---|
| 654 | perihelion = semimaj_axis*(1._dp - eccentricity) |
|---|
| 655 | aphelion = semimaj_axis*(1._dp + eccentricity) |
|---|
| 656 | |
|---|
| 657 | call print_msg('Year (ini|now) = '//real2str(old_year)//' | '//real2str(current_year)) |
|---|
| 658 | call print_msg('Obl. (ini|now) = '//real2str(obl_old)//' | '//real2str(obliquity)) |
|---|
| 659 | call print_msg('Ecc. (ini|now) = '//real2str(ecc_old)//' | '//real2str(eccentricity)) |
|---|
| 660 | call print_msg('Lsp (ini|now) = '//real2str(lsp_old)//' | '//real2str(ls_peri)) |
|---|
| 661 | |
|---|
| 662 | END SUBROUTINE update_orbit |
|---|
| 663 | !======================================================================= |
|---|
| 664 | |
|---|
| 665 | END MODULE orbit |
|---|
