Context Navigation

orbit_param_criterion_mod.F90 @ 2965

Last change on this file since 2965 was 2963, checked in by romain.vande, 3 years ago

Mars PEM :

Adapt PEM to the subslope PCM configuration, it is now fully compatible.

Create a PEM folder in deftank that contains:

run_pem1: a bash file that runs chained simulation of PEM as well as running a parameterizable number of PCM simulation in between.

It also takes care of reshaping XIOS output as well as renaming outputs etc… in the spirit of run_month1.

It is written for Irene machine and the header needs to be adapted for other machines.

run_PEM.def: A text file that shows the possible parameters to choose before a PEM simulation.

It should be included at the end of run.def just like callphys.def

ob_ex_lsp.asc: An ascii file containing the obliquity, eccentricity, ls_peri data from Laskar in Martian year.
README: A txt file explaining the content of the folder

Adapt field_def_physics_mars.xml to consider the case with 7 subslopes in the PCM.
Change context_lmdz_physics.xml to be able to output the file needed by the PEM.

Correct a few other minor bugs.

File size: 7.1 KB

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

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

Download in other formats:

Original Format