Context Navigation

chemthermos.F90 @ 2599

Last change on this file since 2599 was 2042, checked in by emillour, 6 years ago

Mars GCM:
Modifications to use the parametrized photoabsorbtion coefficients;
a first step towards implementing ionospheric chemistry in the new
chemical solver:

change in species indexes in chemthermos.F90, paramfoto_compact.F, hrtherm.F and euvheat.F90
calchim.F90: added a variable in call to photochemistry
photochemistry.F90: added calls to jthermcalc_e107 and phdisrate, with an additionlal flag, jparam (.false. by default). The computed photodissociation coefficents are sent to v_phot, which is used in the chemistry. Thus concentrations computed in chimtogcm are now done over all atmospheric layers.

FGG

File size: 20.2 KB

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

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

Download in other formats:

Original Format