Context Navigation

orbit_param_criterion_mod.F90 @ 2841

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

Mars PEM:
Introduction of the possibility to follow an orbital forcing.
Introduction of new control parameters.
Cleaning of the PEM (removing unused files, add comments and new files)

A file named run_PEM.def can be added to the run.def. It contains the following variables:

_ evol_orbit_pem: Boolean. Do you want to follow an orbital forcing predefined (read in ob_ex_lsp.asc for example)? (default=false)
_ year_bp_ini: Integer. Number of year before present to start the pem run if evol_orbit_pem=.true. , default=0
_ Max_iter_pem: Integer. Maximal number of iteration if none of the stopping criterion is reached and if evol_orbit_pem=.false., default=99999999
_ dt_pem: Integer. Time step of the PEM in year, default=1
_ alpha_criterion: Real. Acceptance rate of sublimating ice surface change, default=0.2
_ soil_pem: Boolean. Do you want to run with subsurface physical processes in the PEM? default=.true.

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

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

Download in other formats:

Original Format