Context Navigation

chemthermos.F90 @ 1268

Last change on this file since 1268 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 19.5 KB

Line
1	SUBROUTINE chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp, &
2	zdens,zpress,zlocal,zenit,ptimestep,zday)
3
4	use tracer_mod, only: nqmx, igcm_co2, igcm_co, igcm_o, igcm_o1d, &
5	igcm_o2, igcm_h, igcm_h2, igcm_oh, igcm_ho2,&
6	igcm_h2o2, igcm_h2o_vap, igcm_o3, igcm_n2, &
7	igcm_n, igcm_no, igcm_no2, igcm_n2d, &
8	igcm_co2plus, igcm_o2plus, igcm_coplus, &
9	igcm_cplus, igcm_nplus, igcm_noplus, &
10	igcm_n2plus, igcm_hplus, igcm_hco2plus, &
11	igcm_elec, igcm_oplus
12	USE comcstfi_h
13	IMPLICIT NONE
14	!=======================================================================
15	! subject:
16	! --------
17	! Computing chemical variations in the thermosphere
18	!
19	! author: MAC July 2003
20	! ------
21	! modifications:
22	! -------------
23	! Ehouarn Sept 2008: added handling of tracers by their names
24	! Francisco Sept 2009: added ionosphere and N chemistry
25	! Francisco Mar 2012: Different chemical schemes possible according to
26	! traceur.def
27	!=======================================================================
28	!
29	! 0. Declarations :
30	! ------------------
31	!
32	#include "callkeys.h"
33	!-----------------------------------------------------------------------
34	! Input/Output
35	! ------------
36	integer :: lswitch,ig,chemthermod,nlayer
37	real :: zycol(nlayer,nqmx)
38	real :: ztemp(nlayer)
39	real :: zdens(nlayer)
40	real :: zpress(nlayer) ! in mbar
41	real :: zlocal(nlayer)
42	real :: zenit
43	real :: ptimestep
44	real :: zday
45	!
46	! Local variables :
47	! -----------------
48	integer :: l
49	integer,save :: nesptherm
50	real,allocatable :: rm(:,:) !number density (cm-3)
51	logical,save :: firstcall=.true.
52
53	! Tracer indexes in the GCM:
54	integer,save :: g_co2=0
55	integer,save :: g_co=0
56	integer,save :: g_o=0
57	integer,save :: g_o1d=0
58	integer,save :: g_o2=0
59	integer,save :: g_o3=0
60	integer,save :: g_h=0
61	integer,save :: g_h2=0
62	integer,save :: g_oh=0
63	integer,save :: g_ho2=0
64	integer,save :: g_h2o2=0
65	integer,save :: g_n2=0
66	integer,save :: g_h2o_vap=0
67	integer,save :: g_n=0
68	integer,save :: g_no=0
69	integer,save :: g_no2=0
70	integer,save :: g_n2d=0
71	integer,save :: g_co2plus=0
72	integer,save :: g_oplus=0
73	integer,save :: g_o2plus=0
74	integer,save :: g_coplus=0
75	integer,save :: g_cplus=0
76	integer,save :: g_nplus=0
77	integer,save :: g_noplus=0
78	integer,save :: g_n2plus=0
79	integer,save :: g_hplus=0
80	integer,save :: g_hco2plus=0
81	integer,save :: g_elec=0
82
83	! Tracer indexes in the thermospheric chemistry:
84	!!! IMPORTANT. These values have to be the same than in
85	!!! jthermcal.F, column.F, paramfoto_compact.F, euvheat.F
86	!!! and hrtherm.F
87	!!! If some value is changed here, the same change has to
88	!!! be made into the other subroutines
89	integer,parameter :: i_co2=1
90	integer,parameter :: i_o2=2
91	integer,parameter :: i_o=3
92	integer,parameter :: i_co=4
93	integer,parameter :: i_h=5
94	integer,parameter :: i_oh=6
95	integer,parameter :: i_ho2=7
96	integer,parameter :: i_h2=8
97	integer,parameter :: i_h2o=9
98	integer,parameter :: i_h2o2=10
99	integer,parameter :: i_o1d=11
100	integer,parameter :: i_o3=12
101	integer,parameter :: i_n2=13
102	integer,parameter :: i_n=14
103	integer,parameter :: i_no=15
104	integer,parameter :: i_n2d=16
105	integer,parameter :: i_no2=17
106	integer,parameter :: i_co2plus=18
107	integer,parameter :: i_oplus=19
108	integer,parameter :: i_o2plus=20
109	integer,parameter :: i_coplus=21
110	integer,parameter :: i_cplus=22
111	integer,parameter :: i_nplus=23
112	integer,parameter :: i_noplus=24
113	integer,parameter :: i_n2plus=25
114	integer,parameter :: i_hplus=26
115	integer,parameter :: i_hco2plus=27
116	integer,parameter :: i_elec=28
117
118
119	! Initializations at first call
120	if (firstcall) then
121	!Number of species to consider
122	nesptherm=0
123	! get the indexes of the tracers we'll need
124	! It is not really necessary to check if the species are present,
125	! as it is already done in calchim.F90. But we include it to double-check
126
127	! Species always included if photochemistry
128	g_co2=igcm_co2
129	if (g_co2.eq.0) then
130	write(,) "chemthermos: Error; no CO2 tracer !!!"
131	write(,) "CO2 must always be present if thermochem=T"
132	stop
133	else
134	nesptherm=nesptherm+1
135	endif
136	g_co=igcm_co
137	if (g_co.eq.0) then
138	write(,) "chemthermos: Error; no CO tracer !!!"
139	write(,) "CO must always be present if thermochem=T"
140	stop
141	else
142	nesptherm=nesptherm+1
143	endif
144	g_o=igcm_o
145	if (g_o.eq.0) then
146	write(,) "chemthermos: Error; no O tracer !!!"
147	write(,) "O must always be present if thermochem=T"
148	stop
149	else
150	nesptherm=nesptherm+1
151	endif
152	g_o1d=igcm_o1d
153	if (g_o1d.eq.0) then
154	write(,) "chemthermos: Error; no O1D tracer !!!"
155	write(,) "O1D must always be present if thermochem=T"
156	stop
157	else
158	nesptherm=nesptherm+1
159	endif
160	g_o2=igcm_o2
161	if (g_o2.eq.0) then
162	write(,) "chemthermos: Error; no O2 tracer !!!"
163	write(,) "O2 must always be present if thermochem=T"
164	stop
165	else
166	nesptherm=nesptherm+1
167	endif
168	g_h=igcm_h
169	if (g_h.eq.0) then
170	write(,) "chemthermos: Error; no H tracer !!!"
171	write(,) "H must always be present if thermochem=T"
172	stop
173	else
174	nesptherm=nesptherm+1
175	endif
176	g_h2=igcm_h2
177	if (g_h2.eq.0) then
178	write(,) "chemthermos: Error; no H2 tracer !!!"
179	write(,) "H2 must always be present if thermochem=T"
180	stop
181	else
182	nesptherm=nesptherm+1
183	endif
184	g_oh=igcm_oh
185	if (g_oh.eq.0) then
186	write(,) "chemthermos: Error; no OH tracer !!!"
187	write(,) "OH must always be present if thermochem=T"
188	stop
189	else
190	nesptherm=nesptherm+1
191	endif
192	g_ho2=igcm_ho2
193	if (g_ho2.eq.0) then
194	write(,) "chemthermos: Error; no HO2 tracer !!!"
195	write(,) "HO2 must always be present if thermochem=T"
196	stop
197	else
198	nesptherm=nesptherm+1
199	endif
200	g_h2o2=igcm_h2o2
201	if (g_h2o2.eq.0) then
202	write(,) "chemthermos: Error; no H2O2 tracer !!!"
203	write(,) "H2O2 must always be present if thermochem=T"
204	stop
205	else
206	nesptherm=nesptherm+1
207	endif
208	g_h2o_vap=igcm_h2o_vap
209	if (g_h2o_vap.eq.0) then
210	write(,) "chemthermos: Error; no water vapor tracer !!!"
211	write(,) "H2O must always be present if thermochem=T"
212	stop
213	else
214	nesptherm=nesptherm+1
215	endif
216	!Verify if O3 chemistry wanted
217	g_o3=igcm_o3
218	if(chemthermod.ge.1) then
219	if (g_o3.eq.0) then
220	write(,) "chemthermos: Error; no O3 tracer !!!"
221	write(,) "O3 must be present if NO is in traceur.def"
222	stop
223	else
224	nesptherm=nesptherm+1
225	endif
226	else
227	if(g_o3 == 0) then
228	write(,) "chemthermos: No O3 chemistry"
229	else
230	write(,) "chemthermos: O3 chemistry included"
231	nesptherm=nesptherm+1
232	chemthermod=1
233	endif
234	endif ! Of if(chemthermod.ge.1)
235
236	! N chemistry
237	if(chemthermod.ge.2) then
238	g_n2=igcm_n2
239	if (g_n2.eq.0) then
240	write(,) "chemthermos: Error; no N2 tracer !!!"
241	write(,) "N2 must be present if NO is in traceur.def"
242	stop
243	else
244	nesptherm=nesptherm+1
245	endif
246	g_n=igcm_n
247	if (g_n.eq.0) then
248	write(,) "chemthermos: Error; no N tracer !!!"
249	write(,) "N must be present if NO is in traceur.def"
250	stop
251	else
252	nesptherm=nesptherm+1
253	endif
254	g_no=igcm_no
255	if (g_no.eq.0) then
256	write(,) "chemthermos: Error; no NO tracer !!!"
257	write(,) "NO must be present if N chemistry wanted"
258	stop
259	else
260	nesptherm=nesptherm+1
261	endif
262	g_no2=igcm_no2
263	if (g_no2.eq.0) then
264	write(,) "chemthermos: Error; no NO2 tracer !!!"
265	write(,) "NO must be present if NO is in traceur.def"
266	stop
267	else
268	nesptherm=nesptherm+1
269	endif
270	g_n2d=igcm_n2d
271	if (g_n2d.eq.0) then
272	write(,) "chemthermos: Error; no N2D tracer !!!"
273	write(,) "NO must be present if NO is in traceur.def"
274	stop
275	else
276	nesptherm=nesptherm+1
277	endif
278	endif ! Of if(chemthermod.ge.2)
279
280	!Ion chemistry
281	if(chemthermod.eq.3) then
282	g_co2plus=igcm_co2plus
283	if (g_co2plus.eq.0 .and. chemthermod.ge.3) then
284	write(,) "chemthermos: Error; no CO2+ tracer !!!"
285	write(,) "CO2+ must be present if chemthermod=3"
286	stop
287	else
288	nesptherm=nesptherm+1
289	endif
290	g_oplus=igcm_oplus
291	if (g_oplus.eq.0 .and. chemthermod.ge.3) then
292	write(,) "chemthermos: Error; no O+ tracer !!!"
293	write(,) "O+ must be present if chemthermod=3"
294	stop
295	else
296	nesptherm=nesptherm+1
297	endif
298	g_o2plus=igcm_o2plus
299	if (g_o2plus.eq.0 .and. chemthermod.ge.3) then
300	write(,) "chemthermos: Error; no O2+ tracer !!!"
301	write(,) "O2+ must be present if chemthermod=3"
302	stop
303	else
304	nesptherm=nesptherm+1
305	endif
306	g_coplus=igcm_coplus
307	if (g_coplus.eq.0) then
308	write(,) "chemthermos: Error; no CO+ tracer !!!"
309	write(,) "CO+ must be present if chemthermod=3"
310	stop
311	else
312	nesptherm=nesptherm+1
313	endif
314	g_cplus=igcm_cplus
315	if (g_cplus.eq.0) then
316	write(,) "chemthermos: Error; no C+ tracer !!!"
317	write(,) "C+ must be present if chemthermod=3"
318	stop
319	else
320	nesptherm=nesptherm+1
321	endif
322	g_nplus=igcm_nplus
323	if (g_nplus.eq.0) then
324	write(,) "chemthermos: Error; no N+ tracer !!!"
325	write(,) "N+ must be present if chemthermod=3"
326	stop
327	else
328	nesptherm=nesptherm+1
329	endif
330	g_noplus=igcm_noplus
331	if (g_noplus.eq.0) then
332	write(,) "chemthermos: Error; no NO+ tracer !!!"
333	write(,) "NO+ must be present if chemthermod=3"
334	stop
335	else
336	nesptherm=nesptherm+1
337	endif
338	g_n2plus=igcm_n2plus
339	if (g_n2plus.eq.0) then
340	write(,) "chemthermos: Error; no N2+ tracer !!!"
341	write(,) "N2+ must be present if chemthermod=3"
342	stop
343	else
344	nesptherm=nesptherm+1
345	endif
346	g_hplus=igcm_hplus
347	if (g_hplus.eq.0) then
348	write(,) "chemthermos: Error; no H+ tracer !!!"
349	write(,) "H+ must be present if chemthermod=3"
350	stop
351	else
352	nesptherm=nesptherm+1
353	endif
354	g_hco2plus=igcm_hco2plus
355	if (g_hco2plus.eq.0 .and. chemthermod.ge.3) then
356	write(,) "chemthermos: Error; no HCO2+ tracer !!!"
357	write(,) "HCO2+ must be present if chemthermod=3"
358	stop
359	else
360	nesptherm=nesptherm+1
361	endif
362	g_elec=igcm_elec
363	if (g_elec.eq.0) then
364	write(,) "chemthermos: Error; no e- tracer !!!"
365	write(,) "e- must be present if chemthermod=3"
366	stop
367	else
368	nesptherm=nesptherm+1
369	endif
370	endif
371
372	!Check if number of species is as expected
373	select case(chemthermod)
374	case(0)
375	if(nesptherm.ne.11) then
376	write(,)"chemthermos: Error:"
377	write(,)"Number of tracers not correct"
378	write(,)"There are",nesptherm," tracers"
379	write(,)"There should be 11 tracers"
380	stop
381	else
382	write(,)"chemthermos:",nesptherm," species"
383	endif
384	case(1)
385	if(nesptherm.ne.12) then
386	write(,)"chemthermos: Error:"
387	write(,)"Number of tracers not correct"
388	write(,)"There are",nesptherm," tracers"
389	write(,)"There should be 12 tracers"
390	stop
391	else
392	write(,)"chemthermos:",nesptherm," species"
393	endif
394	case(2)
395	if(nesptherm.ne.17) then
396	write(,)"chemthermos: Error:"
397	write(,)"Number of tracers not correct"
398	write(,)"There are",nesptherm," tracers"
399	write(,)"There should be 17 tracers"
400	stop
401	else
402	write(,)"chemthermos:",nesptherm," species"
403	endif
404	case(3)
405	if(nesptherm.ne.28) then
406	write(,)"chemthermos: Error:"
407	write(,)"Number of tracers not correct"
408	write(,)"There are",nesptherm," tracers"
409	write(,)"There should be 28 tracers"
410	stop
411	else
412	write(,)"chemthermos:",nesptherm," species"
413	endif
414	end select
415	firstcall= .false.
416	write(,)'chemthermos: chemthermod=',chemthermod
417	endif ! of if (firstcall)
418
419	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
420
421	!Allocate density vector
422	allocate(rm(nlayer,nesptherm))
423
424	do l=1,nlayer
425	rm(l,i_co2) = max(zycol(l,g_co2)*zdens(l),1.e-30)
426	rm(l,i_co) = max(zycol(l,g_co)*zdens(l),1.e-30)
427	rm(l,i_o) = max(zycol(l,g_o)*zdens(l),1.e-30)
428	rm(l,i_o1d) = max(zycol(l,g_o1d)*zdens(l),1.e-30)
429	rm(l,i_o2) = max(zycol(l,g_o2)*zdens(l),1.e-30)
430	rm(l,i_h) = max(zycol(l,g_h)*zdens(l),1.e-30)
431	rm(l,i_h2) = max(zycol(l,g_h2)*zdens(l),1.e-30)
432	rm(l,i_oh) = max(zycol(l,g_oh)*zdens(l),1.e-30)
433	rm(l,i_ho2) = max(zycol(l,g_ho2)*zdens(l),1.e-30)
434	rm(l,i_h2o2) = max(zycol(l,g_h2o2)*zdens(l),1.e-30)
435	rm(l,i_h2o) = max(zycol(l,g_h2o_vap)*zdens(l),1.e-30)
436	if(chemthermod.ge.1) &
437	rm(l,i_o3) = max(zycol(l,g_o3)*zdens(l),1.e-30)
438	if(chemthermod.ge.2) then
439	rm(l,i_n2) = max(zycol(l,g_n2)*zdens(l),1.e-30)
440	rm(l,i_n) = max(zycol(l,g_n)*zdens(l),1.e-30)
441	rm(l,i_no) = max(zycol(l,g_no)*zdens(l),1.e-30)
442	rm(l,i_co) = max(zycol(l,g_co)*zdens(l),1.e-30)
443	rm(l,i_n2d) = max(zycol(l,g_n2d)*zdens(l),1.e-30)
444	rm(l,i_no2) = max(zycol(l,g_no2)*zdens(l),1.e-30)
445	endif
446	if(chemthermod.eq.3) then
447	rm(l,i_co2plus) = max(zycol(l,g_co2plus)*zdens(l),1.e-30)
448	rm(l,i_oplus) = max(zycol(l,g_oplus)*zdens(l),1.e-30)
449	rm(l,i_o2plus) = max(zycol(l,g_o2plus)*zdens(l),1.e-30)
450	rm(l,i_coplus) = max(zycol(l,g_coplus)*zdens(l),1.e-30)
451	rm(l,i_cplus) = max(zycol(l,g_cplus)*zdens(l),1.e-30)
452	rm(l,i_nplus) = max(zycol(l,g_nplus)*zdens(l),1.e-30)
453	rm(l,i_noplus) = max(zycol(l,g_noplus)*zdens(l),1.e-30)
454	rm(l,i_n2plus) = max(zycol(l,g_n2plus)*zdens(l),1.e-30)
455	rm(l,i_hplus) = max(zycol(l,g_hplus)*zdens(l),1.e-30)
456	rm(l,i_hco2plus)= max(zycol(l,g_hco2plus)*zdens(l),1.e-30)
457	rm(l,i_elec) = max(zycol(l,g_elec)*zdens(l),1.e-30)
458	endif
459	if(chemthermod.gt.3.or.chemthermod.lt.0) then
460	write(,)'chemthermos: bad value for chemthermod. Stop'
461	stop
462	endif
463	enddo
464
465	!Solar flux calculation
466
467	!Photoabsorption coefficients
468	call flujo(solarcondate)
469	if(solvarmod.eq.0) then
470	call jthermcalc(ig,nlayer,chemthermod,rm,nesptherm,ztemp,zlocal,zenit)
471	else if(solvarmod.eq.1) then
472	call jthermcalc_e107(ig,nlayer,chemthermod,rm,nesptherm,ztemp,zlocal,zenit,zday)
473	endif
474
475
476	!Chemistry
477	call paramfoto_compact &
478	(ig,nlayer,chemthermod,lswitch,ztemp,ptimestep,zenit,zlocal,rm,nesptherm)
479
480	!Concentrations back to GCM
481	do l=lswitch,nlayer
482	zycol(l,g_co2) = max(rm(l,i_co2) / zdens(l) , 1.e-30)
483	zycol(l,g_co) = max(rm(l,i_co) / zdens(l) , 1.e-30)
484	zycol(l,g_o2) = max(rm(l,i_o2) / zdens(l) , 1.e-30)
485	zycol(l,g_h2) = max(rm(l,i_h2) / zdens(l) , 1.e-30)
486	zycol(l,g_h) = max(rm(l,i_h) / zdens(l) , 1.e-30)
487	zycol(l,g_oh) = max(rm(l,i_oh) / zdens(l) , 1.e-30)
488	zycol(l,g_ho2) = max(rm(l,i_ho2) / zdens(l) , 1.e-30)
489	zycol(l,g_h2o_vap) = max(rm(l,i_h2o) / zdens(l) , 1.e-30)
490	zycol(l,g_h2o2) = max(rm(l,i_h2o2) / zdens(l) , 1.e-30)
491	zycol(l,g_o1d) = max(rm(l,i_o1d) / zdens(l) , 1.e-30)
492	zycol(l,g_o) = max(rm(l,i_o) / zdens(l) , 1.e-30)
493	if(chemthermod.ge.1) &
494	zycol(l,g_o3) = max(rm(l,i_o3) / zdens(l) , 1.e-30)
495	if(chemthermod.ge.2) then
496	zycol(l,g_n2) = max(rm(l,i_n2) / zdens(l) , 1.e-30)
497	zycol(l,g_n) = max(rm(l,i_n) / zdens(l) , 1.e-30)
498	zycol(l,g_no) = max(rm(l,i_no) / zdens(l) , 1.e-30)
499	zycol(l,g_no2) = max(rm(l,i_no2) / zdens(l) , 1.e-30)
500	zycol(l,g_n2d) = max(rm(l,i_n2d) / zdens(l) , 1.e-30)
501	endif
502	if(chemthermod.ge.3) then
503	zycol(l,g_co2plus) = max(rm(l,i_co2plus) / zdens(l) , 1.e-30)
504	zycol(l,g_oplus) = max(rm(l,i_oplus) / zdens(l) , 1.e-30)
505	zycol(l,g_o2plus) = max(rm(l,i_o2plus) / zdens(l) , 1.e-30)
506	zycol(l,g_coplus) = max(rm(l,i_coplus) / zdens(l) , 1.e-30)
507	zycol(l,g_cplus) = max(rm(l,i_cplus) / zdens(l) , 1.e-30)
508	zycol(l,g_nplus) = max(rm(l,i_nplus) / zdens(l) , 1.e-30)
509	zycol(l,g_noplus) = max(rm(l,i_noplus) / zdens(l) , 1.e-30)
510	zycol(l,g_n2plus) = max(rm(l,i_n2plus) / zdens(l) , 1.e-30)
511	zycol(l,g_hplus) = max(rm(l,i_hplus) / zdens(l) , 1.e-30)
512	zycol(l,g_hco2plus)= max(rm(l,i_hco2plus)/ zdens(l) , 1.e-30)
513	zycol(l,g_elec) = max(rm(l,i_elec) / zdens(l) , 1.e-30)
514	endif
515	enddo !nlayer
516
517	!Deallocations
518	deallocate(rm)
519
520	return
521	end
522
523
524
525

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

Download in other formats:

Original Format