Context Navigation

orbit_param_criterion_mod.F90 @ 3019

Last change on this file since 3019 was 2986, checked in by romain.vande, 2 years ago
MARS PEM: Small fix to adapt pem to the r2981 write_profile modification by evos. RV
File size: 8.1 KB

Line
1	MODULE orbit_param_criterion_mod
2
3	!=======================================================================
4	!
5	! Purpose: Compute the maximum nimber of iteration the PEM can do based
6	! on the stopping criterion given by the orbital parameters
7	!
8	! Author: RV
9	!=======================================================================
10
11	IMPLICIT NONE
12
13	CONTAINS
14
15	SUBROUTINE orbit_param_criterion(year_iter_max)
16	#ifdef CPP_IOIPSL
17	use IOIPSL, only: getin
18	#else
19	! if not using IOIPSL, we still need to use (a local version of) getin
20	use ioipsl_getincom, only: getin
21	#endif
22
23	USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
24	#ifndef CPP_STD
25	USE planete_h, ONLY: e_elips, obliquit, timeperi
26	#else
27	use planete_mod, only: e_elips, obliquit, timeperi
28	#endif
29	USE comconst_mod, ONLY: pi
30	USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
31	ini_lask_param_mod,last_ilask
32
33	IMPLICIT NONE
34
35	!--------------------------------------------------------
36	! Input Variables
37	!--------------------------------------------------------
38
39	!--------------------------------------------------------
40	! Output Variables
41	!--------------------------------------------------------
42
43	integer,intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
44
45	!--------------------------------------------------------
46	! Local variables
47	!--------------------------------------------------------
48
49	real :: Year ! Year of the simulation
50	real :: timeperi_ls ! Year of the simulation
51	integer nlask,ilask!,last_ilask ! Loop variable
52	parameter (nlask = 20001) ! Number of line in Laskar file
53
54	real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible
55
56	real max_obl,max_ex,max_lsp ! Maximum value of orbit param given the acceptable percentage
57	real min_obl,min_ex,min_lsp ! Maximum value of orbit param given the acceptable percentage
58
59	real max_obl_iter,max_ex_iter,max_lsp_iter ! Maximum year iteration before reaching an unacceptable value
60
61	logical obl_not_found, ex_not_found,lsp_not_found ! Loop variable (first call)
62	logical max_lsp_modulo,min_lsp_modulo ! Variable to check if we are close to the 360 degree modulo for lsp parameter
63
64	real xa,xb,ya,yb,yc
65
66	! **********************************************************************
67	! 0. Initializations
68	! **********************************************************************
69
70	Year=year_bp_ini+year_PEM ! We compute the current year
71	timeperi_ls=timeperi*360/2/pi ! We convert in degree
72
73	call ini_lask_param_mod(nlask) ! Allocation of variables
74
75	print *, "orbit_param_criterion, Year in pem.def=", year_bp_ini
76	print *, "orbit_param_criterion, Year in the startpem.nc =", year_PEM
77	print *, "orbit_param_criterion, Current year=", Year
78	print *, "orbit_param_criterion, Current obl=", obliquit
79	print *, "orbit_param_criterion, Current ex=", e_elips
80	print *, "orbit_param_criterion, Current lsp=", timeperi_ls
81
82	! We read the file
83
84	open(73,file='ob_ex_lsp.asc')
85	do ilask=1,nlask
86	read(73,*) yearlask(ilask),oblask(ilask), &
87	exlask(ilask),lsplask(ilask)
88	yearlask(ilask)=yearlask(ilask)*1000
89
90	if(yearlask(ilask).GT.Year) then
91	last_ilask=ilask+1
92	endif
93	end do
94	close(73)
95
96	print *, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask
97
98	! 5% max change case
99
100	max_change_obl=0.05
101	max_change_ex=0.05
102	max_change_lsp=0.05
103
104	max_obl=obliquit*(1+max_change_obl)
105	min_obl=obliquit*(1-max_change_obl)
106
107	max_ex=e_elips*(1+max_change_ex)
108	min_ex=e_elips*(1-max_change_ex)
109
110	max_lsp=timeperi_ls*(1+max_change_lsp)
111	min_lsp=timeperi_ls*(1-max_change_lsp)
112
113	!End of 5% max change case
114
115	!Constant max change case
116
117	max_change_obl=0.5
118	max_change_ex=0.1
119	max_change_lsp=40.
120
121	CALL getin('max_change_obl', max_change_obl)
122
123	CALL getin('max_change_ex', max_change_ex)
124
125	CALL getin('max_change_lsp', max_change_lsp)
126
127	max_obl=obliquit+max_change_obl
128	min_obl=obliquit-max_change_obl
129
130	max_ex=e_elips+max_change_ex
131	min_ex=e_elips-max_change_ex
132
133	max_lsp=timeperi_ls+max_change_lsp
134	min_lsp=timeperi_ls-max_change_lsp
135
136	max_lsp_modulo=.false.
137	min_lsp_modulo=.false.
138	if(max_lsp.gt.360) max_lsp_modulo=.true.
139	if(min_lsp.lt.0) min_lsp_modulo=.true.
140
141	!End Constant max change case
142
143	! If we do not want some orb parameter to change, they should not be a stopping criterion,
144	! So the number of iteration corresponding is set to maximum
145	if(.not.var_obl) then
146	obl_not_found=.FALSE.
147	else
148	obl_not_found=.TRUE.
149	endif
150	if(.not.var_ex) then
151	ex_not_found=.FALSE.
152	else
153	ex_not_found=.TRUE.
154	endif
155	if(.not.var_lsp) then
156	lsp_not_found=.FALSE.
157	else
158	lsp_not_found=.TRUE.
159	endif
160
161	max_obl_iter=999999
162	max_ex_iter =999999
163	max_lsp_iter=999999
164
165	!--------------------------------
166
167	yc=max_obl
168	yc=min_obl
169
170	!--------------------------------
171
172	do ilask=last_ilask,2,-1
173	xa=yearlask(ilask)
174	xb=yearlask(ilask-1)
175	! Linear interpolation to find the maximum number of iteration
176	if((oblask(ilask-1).GT.max_obl).and. obl_not_found ) then
177	ya=oblask(ilask)
178	yb=oblask(ilask-1)
179	yc=max_obl
180	max_obl_iter=xa(yc-yb)/(ya-yb)+xb(ya-yc)/(ya-yb)-Year
181	obl_not_found=.FALSE.
182	elseif((oblask(ilask-1).LT.min_obl).and. obl_not_found ) then
183	ya=oblask(ilask)
184	yb=oblask(ilask-1)
185	yc=min_obl
186	max_obl_iter=xa(yc-yb)/(ya-yb)+xb(ya-yc)/(ya-yb)-Year
187	obl_not_found=.FALSE.
188	endif
189	if((exlask(ilask-1).GT.max_ex).and. ex_not_found ) then
190	ya=exlask(ilask)
191	yb=exlask(ilask-1)
192	yc=max_ex
193	max_ex_iter=xa(yc-yb)/(ya-yb)+xb(ya-yc)/(ya-yb)-Year
194	ex_not_found=.FALSE.
195	elseif((exlask(ilask-1).LT.min_ex ).and. ex_not_found ) then
196	ya=exlask(ilask)
197	yb=exlask(ilask-1)
198	yc=min_ex
199	max_ex_iter=xa(yc-yb)/(ya-yb)+xb(ya-yc)/(ya-yb)-Year
200	ex_not_found=.FALSE.
201	endif
202	if(max_lsp_modulo .and. lsplask(ilask).lt.180) lsplask(ilask)=lsplask(ilask)+360.
203	if(min_lsp_modulo .and. lsplask(ilask).gt.180) lsplask(ilask)=lsplask(ilask)-360.
204	if((lsplask(ilask-1).GT.max_lsp).and. lsp_not_found ) then
205	ya=lsplask(ilask)
206	yb=lsplask(ilask-1)
207	yc=max_lsp
208	max_ex_iter=xa(yc-yb)/(ya-yb)+xb(ya-yc)/(ya-yb)-Year
209	lsp_not_found=.FALSE.
210	elseif((lsplask(ilask-1).LT.min_lsp ).and. lsp_not_found ) then
211	ya=lsplask(ilask)
212	yb=lsplask(ilask-1)
213	yc=min_lsp
214	max_ex_iter=xa(yc-yb)/(ya-yb)+xb(ya-yc)/(ya-yb)-Year
215	lsp_not_found=.FALSE.
216	endif
217	enddo
218
219	print *, "Maximum obliquity accepted=", max_obl
220	print *, "Minimum obliquity accepted=", min_obl
221	print *, "Maximum number of iteration for the obl. parameter=", max_obl_iter
222
223	print *, "Maximum excentricity accepted=", max_ex
224	print *, "Minimum excentricity accepted=", min_ex
225	print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter
226
227	print *, "Maximum lsp accepted=", max_lsp
228	print *, "Minimum lsp accepted=", min_lsp
229	print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter
230
231	year_iter_max=min(max_obl_iter,max_ex_iter,max_lsp_iter)
232
233	print *, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max
234
235	END SUBROUTINE orbit_param_criterion
236
237	!********************************************************************************
238
239	END MODULE orbit_param_criterion_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format