source: trunk/LMDZ.COMMON/libf/evolution/orbit.F90 @ 4117

Last change on this file since 4117 was 4117, checked in by jbclement, 2 days ago

PEM:

  • Remove 'd_h2oice' and 'd_co2ice' arguments from adsorption and soil thermal inertia routines (not needed and possible bugs).
  • Move 'read_callphys' call from 'read_start1D' to 'read_rundef' so that 'CO2cond_ps' is initialized no matter what.
  • Explicit initialization of 'delta_h2o_ads'/'delta_co2_ads' before 'read_startpem'.
  • Add time index argument to 'put_var_nc' calls to write "restart.nc".
  • Increase temporary string buffer size for 'real2str_qp'.
  • Fixed some typos.

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