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 | logical :: crossed_modulo_d, crossed_modulo_i ! Flag variables for modulo of Lsp |
---|
60 | |
---|
61 | ! ********************************************************************** |
---|
62 | ! 0. Initializations |
---|
63 | ! ********************************************************************** |
---|
64 | Year = year_bp_ini + i_myear ! We compute the current year |
---|
65 | Lsp = lsperi*180./pi ! We convert in degree |
---|
66 | |
---|
67 | call ini_lask_param_mod ! Allocation of variables |
---|
68 | |
---|
69 | write(*,*) 'orbit_param_criterion, Current year =', Year |
---|
70 | write(*,*) 'orbit_param_criterion, Current obl. =', obliquit |
---|
71 | write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips |
---|
72 | write(*,*) 'orbit_param_criterion, Current Lsp =', Lsp |
---|
73 | |
---|
74 | ! We read the file |
---|
75 | open(73,file = 'obl_ecc_lsp.asc') |
---|
76 | last_ilask = 0 |
---|
77 | do ilask = 1,size(yearlask,1) |
---|
78 | read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask) |
---|
79 | yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years |
---|
80 | if (yearlask(ilask) > Year) last_ilask = ilask + 1 |
---|
81 | enddo |
---|
82 | close(73) |
---|
83 | if (last_ilask == 0 .or. last_ilask == size(yearlask,1) + 1) then |
---|
84 | error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".' |
---|
85 | else |
---|
86 | write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', last_ilask |
---|
87 | endif |
---|
88 | |
---|
89 | ! Constant max change case |
---|
90 | max_change_obl = 0.6 |
---|
91 | max_change_ecc = 2.8e-3 |
---|
92 | max_change_lsp = 18. |
---|
93 | |
---|
94 | call getin('max_change_obl',max_change_obl) |
---|
95 | max_obl = obliquit + max_change_obl |
---|
96 | min_obl = obliquit - max_change_obl |
---|
97 | write(*,*) 'Maximum obliquity accepted =', max_obl |
---|
98 | write(*,*) 'Minimum obliquity accepted =', min_obl |
---|
99 | |
---|
100 | call getin('max_change_ecc',max_change_ecc) |
---|
101 | max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ecc < 1. |
---|
102 | min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 0. |
---|
103 | write(*,*) 'Maximum eccentricity accepted =', max_ecc |
---|
104 | write(*,*) 'Minimum eccentricity accepted =', min_ecc |
---|
105 | |
---|
106 | call getin('max_change_lsp',max_change_lsp) |
---|
107 | max_lsp = Lsp + max_change_lsp |
---|
108 | min_lsp = Lsp - max_change_lsp |
---|
109 | write(*,*) 'Maximum Lsp accepted =', max_lsp |
---|
110 | write(*,*) 'Minimum Lsp accepted =', min_lsp |
---|
111 | ! End Constant max change case |
---|
112 | |
---|
113 | ! If we do not want some orb parameter to change, they should not be a stopping criterion, |
---|
114 | ! So the number of iteration corresponding is set to maximum |
---|
115 | max_obl_iter = 1000000. |
---|
116 | max_ecc_iter = 1000000. |
---|
117 | max_lsp_iter = 1000000. |
---|
118 | |
---|
119 | ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters |
---|
120 | found_obl = .false. |
---|
121 | found_ecc = .false. |
---|
122 | found_lsp = .false. |
---|
123 | if (.not. var_obl) found_obl = .true. |
---|
124 | if (.not. var_ecc) found_ecc = .true. |
---|
125 | if (.not. var_lsp) found_lsp = .true. |
---|
126 | crossed_modulo_d = .false. |
---|
127 | crossed_modulo_i = .false. |
---|
128 | |
---|
129 | do ilask = last_ilask,2,-1 |
---|
130 | xa = yearlask(ilask) |
---|
131 | xb = yearlask(ilask - 1) |
---|
132 | |
---|
133 | ! Obliquity |
---|
134 | if (var_obl .and. .not. found_obl) then |
---|
135 | ya = obllask(ilask) |
---|
136 | yb = obllask(ilask - 1) |
---|
137 | if (yb < min_obl) then ! The minimum accepted is overtaken |
---|
138 | max_obl_iter = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year |
---|
139 | found_obl = .true. |
---|
140 | else if (max_obl < yb) then ! The maximum accepted is overtaken |
---|
141 | max_obl_iter = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year |
---|
142 | found_obl = .true. |
---|
143 | endif |
---|
144 | endif |
---|
145 | |
---|
146 | ! Eccentricity |
---|
147 | if (var_ecc .and. .not. found_ecc) then |
---|
148 | ya = ecclask(ilask) |
---|
149 | yb = ecclask(ilask - 1) |
---|
150 | if (yb < min_ecc) then ! The minimum accepted is overtaken |
---|
151 | max_ecc_iter = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year |
---|
152 | found_ecc = .true. |
---|
153 | else if (max_ecc < yb) then ! The maximum accepted is overtaken |
---|
154 | max_ecc_iter = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year |
---|
155 | found_ecc = .true. |
---|
156 | endif |
---|
157 | endif |
---|
158 | |
---|
159 | ! Lsp |
---|
160 | if (var_lsp .and. .not. found_lsp) then |
---|
161 | ya = lsplask(ilask) |
---|
162 | yb = lsplask(ilask - 1) |
---|
163 | if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval |
---|
164 | if (ya < yb) then ! Lsp should be decreasing |
---|
165 | if (.not. crossed_modulo_i) then |
---|
166 | yb = yb - 360. |
---|
167 | if (min_lsp < 0.) crossed_modulo_d = .true. |
---|
168 | else |
---|
169 | crossed_modulo_i = .false. |
---|
170 | endif |
---|
171 | else ! Lsp should be increasing |
---|
172 | if (.not. crossed_modulo_d) then |
---|
173 | yb = yb + 360. |
---|
174 | if (max_lsp > 360.) crossed_modulo_i = .true. |
---|
175 | else |
---|
176 | crossed_modulo_d = .false. |
---|
177 | endif |
---|
178 | endif |
---|
179 | else |
---|
180 | if (crossed_modulo_d) then ! If decreasing Lsp "crossed" the modulo but min_lsp has not been met yet |
---|
181 | ya = ya - 360. |
---|
182 | yb = yb - 360. |
---|
183 | else if (crossed_modulo_i) then ! If increasing Lsp "crossed" the modulo but max_lsp has not been met yet |
---|
184 | ya = ya + 360. |
---|
185 | yb = yb + 360. |
---|
186 | endif |
---|
187 | endif |
---|
188 | if (yb < min_lsp) then ! The minimum accepted is overtaken |
---|
189 | max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year |
---|
190 | found_lsp = .true. |
---|
191 | else if (max_lsp < yb) then ! The maximum accepted is overtaken |
---|
192 | max_lsp_iter = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year |
---|
193 | found_lsp = .true. |
---|
194 | endif |
---|
195 | endif |
---|
196 | |
---|
197 | if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found |
---|
198 | enddo |
---|
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 | |
---|