Context Navigation

chemthermos.F90 @ 674

Last change on this file since 674 was 635, checked in by emillour, 13 years ago

Mars GCM: Update of the chemistry package, including:

93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
computation of photoabsorption using raytracing
improved time stepping in the photochemistry
updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
Upadted surfacearea.F
Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM

File size: 18.9 KB

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

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

Download in other formats:

Original Format