Context Navigation

chemthermos.F90 @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

File size: 19.7 KB

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

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

Download in other formats:

Original Format