Context Navigation

calchim_asis.F90 @ 1888

Last change on this file since 1888 was 1888, checked in by emillour, 7 years ago

Mars GCM:
Updated the calculation of the dissociation and ionization
branching ratios, using the values in the Schunk and Nagy book.
New datafiles *branchingratio_schunkandnagy2000_param.dat must be loaded
for the "EUVDAT" subdirectory of the standard "datadir" directory.
Main changes are:

param_v4_h.F90 -> New definition for the O2 ionization branching ratio
param_read_e107.F -> Read the new files containing the S&N branching ratios
paramfoto_compact.F -> Mainly cleaning of the code and the comments.

Also correction of a bug affecting the calculation of CO losses

chemthermos.F90 -> Small modification to add the possibility of including

NO and O2 nightglow rates to the outputs

calchim.F90 and calchim_asis.F90 -> account for change in arguments in

calls to chemthermos
FGG

File size: 27.8 KB

Line
1	subroutine calchim_asis(ngrid,nlayer,nq, &
2	ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, &
3	zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud, &
4	dqscloud,tauref,co2ice, &
5	pu,pdu,pv,pdv,surfdust,surfice)
6
7	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
8	igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2, &
9	igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap, &
10	igcm_no, igcm_n, igcm_no2, igcm_n2d, &
11	igcm_o2plus, igcm_co2plus, igcm_oplus, &
12	igcm_coplus, igcm_cplus, igcm_nplus, &
13	igcm_noplus, igcm_n2plus, igcm_hplus, &
14	igcm_hco2plus, igcm_elec, mmol
15
16	use conc_mod, only: mmean ! mean molecular mass of the atmosphere
17	USE comcstfi_h
18
19	implicit none
20
21	!=======================================================================
22	!
23	! subject:
24	! --------
25	!
26	! Prepare the call for the photochemical module, and send back the
27	! tendencies from photochemistry in the chemical species mass mixing ratios
28	!
29	! Author: Sebastien Lebonnois (08/11/2002)
30	! -------
31	! update 12/06/2003 for water ice clouds and compatibility with dust
32	! update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
33	! update 03/05/2005 cosmetic changes (Franck Lefevre)
34	! update sept. 2008 identify tracers by their names (Ehouarn Millour)
35	! update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
36	! update 16/03/2012 optimization (Franck Lefevre)
37	! update 11/12/2013 optimization (Franck Lefevre)
38	!
39	! Arguments:
40	! ----------
41	!
42	! Input:
43	!
44	! ptimestep timestep (s)
45	! pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
46	! pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa)
47	! pt(ngrid,nlayer) Temperature (K)
48	! pdt(ngrid,nlayer) Temperature tendency (K)
49	! pu(ngrid,nlayer) u component of the wind (ms-1)
50	! pdu(ngrid,nlayer) u component tendency (K)
51	! pv(ngrid,nlayer) v component of the wind (ms-1)
52	! pdv(ngrid,nlayer) v component tendency (K)
53	! dist_sol distance of the sun (AU)
54	! mu0(ngrid) cos of solar zenith angle (=1 when sun at zenith)
55	! pq(ngrid,nlayer,nq) Advected fields, ie chemical species here
56	! pdq(ngrid,nlayer,nq) Previous tendencies on pq
57	! tauref(ngrid) Optical depth at 7 hPa
58	! co2ice(ngrid) co2 ice surface layer (kg.m-2)
59	! surfdust(ngrid,nlayer) dust surface area (m2/m3)
60	! surfice(ngrid,nlayer) ice surface area (m2/m3)
61	!
62	! Output:
63	!
64	! dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
65	! dqschim(ngrid,nq) ! tendencies on qsurf
66	!
67	!=======================================================================
68
69	#include "chimiedata.h"
70	#include "callkeys.h"
71
72	! input:
73
74	integer,intent(in) :: ngrid ! number of atmospheric columns
75	integer,intent(in) :: nlayer ! number of atmospheric layers
76	integer,intent(in) :: nq ! number of tracers
77	real :: ptimestep
78	real :: pplay(ngrid,nlayer) ! pressure at the middle of the layers
79	real :: zzlay(ngrid,nlayer) ! pressure at the middle of the layers
80	real :: pplev(ngrid,nlayer+1) ! intermediate pressure levels
81	real :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
82	real :: pt(ngrid,nlayer) ! temperature
83	real :: pdt(ngrid,nlayer) ! temperature tendency
84	real :: pu(ngrid,nlayer) ! u component of the wind (m.s-1)
85	real :: pdu(ngrid,nlayer) ! u component tendency
86	real :: pv(ngrid,nlayer) ! v component of the wind (m.s-1)
87	real :: pdv(ngrid,nlayer) ! v component tendency
88	real :: dist_sol ! distance of the sun (AU)
89	real :: mu0(ngrid) ! cos of solar zenith angle (=1 when sun at zenith)
90	real :: pq(ngrid,nlayer,nq) ! tracers mass mixing ratio
91	real :: pdq(ngrid,nlayer,nq) ! previous tendencies
92	real :: zday ! date (time since Ls=0, in martian days)
93	real :: tauref(ngrid) ! optical depth at 7 hPa
94	real :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
95	real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
96	real :: surfice(ngrid,nlayer) ! ice surface area (m2/m3)
97
98	! output:
99
100	real :: dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
101	real :: dqschim(ngrid,nq) ! tendencies on qsurf
102	real :: dqcloud(ngrid,nlayer,nq) ! tendencies on pq due to condensation
103	real :: dqscloud(ngrid,nq) ! tendencies on qsurf
104
105	! local variables:
106
107	integer,save :: nbq ! number of tracers used in the chemistry
108	integer,allocatable,save :: niq(:) ! array storing the indexes of the tracers
109	integer :: iloc(1) ! index of major species
110	integer :: ig,l,i,iq,iqmax
111	integer :: foundswitch, lswitch
112	integer,save :: chemthermod
113
114	integer,save :: i_co2 = 0
115	integer,save :: i_co = 0
116	integer,save :: i_o = 0
117	integer,save :: i_o1d = 0
118	integer,save :: i_o2 = 0
119	integer,save :: i_o3 = 0
120	integer,save :: i_h = 0
121	integer,save :: i_h2 = 0
122	integer,save :: i_oh = 0
123	integer,save :: i_ho2 = 0
124	integer,save :: i_h2o2 = 0
125	integer,save :: i_ch4 = 0
126	integer,save :: i_n2 = 0
127	integer,save :: i_h2o = 0
128	integer,save :: i_n = 0
129	integer,save :: i_no = 0
130	integer,save :: i_no2 = 0
131	integer,save :: i_n2d = 0
132	integer,save :: i_co2plus=0
133	integer,save :: i_oplus=0
134	integer,save :: i_o2plus=0
135	integer,save :: i_coplus=0
136	integer,save :: i_cplus=0
137	integer,save :: i_nplus=0
138	integer,save :: i_noplus=0
139	integer,save :: i_n2plus=0
140	integer,save :: i_hplus=0
141	integer,save :: i_hco2plus=0
142	integer,save :: i_elec=0
143
144	integer :: ig_vl1
145
146	real :: latvl1, lonvl1
147	real :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry
148	! new mole fraction after
149	real :: zt(ngrid,nlayer) ! temperature
150	real :: zu(ngrid,nlayer) ! u component of the wind
151	real :: zv(ngrid,nlayer) ! v component of the wind
152	real :: taucol ! optical depth at 7 hPa
153
154	logical,save :: firstcall = .true.
155	logical,save :: depos = .false. ! switch for dry deposition
156
157	! for each column of atmosphere:
158
159	real :: zpress(nlayer) ! Pressure (mbar)
160	real :: zdens(nlayer) ! Density (cm-3)
161	real :: ztemp(nlayer) ! Temperature (K)
162	real :: zlocal(nlayer) ! Altitude (km)
163	real :: zycol(nlayer,nq) ! Composition (mole fractions)
164	real :: zmmean(nlayer) ! Mean molecular mass (g.mole-1)
165	real :: szacol ! Solar zenith angle
166	real :: surfice1d(nlayer) ! Ice surface area (cm2/cm3)
167	real :: surfdust1d(nlayer) ! Dust surface area (cm2/cm3)
168	real :: jo3(nlayer) ! Photodissociation rate O3->O1D (s-1)
169	real :: em_no(nlayer) ! NO nightglow emission rate
170	real :: em_o2(nlayer) ! O2 nightglow emission rate
171
172	integer :: iter(nlayer) ! Number of chemical iterations
173	! within one physical timestep
174
175	! for output:
176
177	logical :: output ! to issue calls to writediagfi and stats
178	parameter (output = .true.)
179	real :: jo3_3d(ngrid,nlayer) ! Photodissociation rate O3->O1D (s-1)
180	real :: emission_no(ngrid,nlayer) !NO emission rate
181	real :: emission_o2(ngrid,nlayer) !O2 emission rate
182	real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
183	! within one physical timestep
184
185	!=======================================================================
186	! initialization of the chemistry (first call only)
187	!=======================================================================
188
189	if (firstcall) then
190
191	if (photochem) then
192	print*,'calchim: Read photolysis lookup table'
193	call read_phototable
194	end if
195	! find index of chemical tracers to use
196	allocate(niq(nq))
197	! Listed here are all tracers that can go into photochemistry
198	nbq = 0 ! to count number of tracers
199	! Species ALWAYS present if photochem=.T. or thermochem=.T.
200	i_co2 = igcm_co2
201	if (i_co2 == 0) then
202	write(,) "calchim: Error; no CO2 tracer !!!"
203	stop
204	else
205	nbq = nbq + 1
206	niq(nbq) = i_co2
207	end if
208	i_co = igcm_co
209	if (i_co == 0) then
210	write(,) "calchim: Error; no CO tracer !!!"
211	stop
212	else
213	nbq = nbq + 1
214	niq(nbq) = i_co
215	end if
216	i_o = igcm_o
217	if (i_o == 0) then
218	write(,) "calchim: Error; no O tracer !!!"
219	stop
220	else
221	nbq = nbq + 1
222	niq(nbq) = i_o
223	end if
224	i_o1d = igcm_o1d
225	if (i_o1d == 0) then
226	write(,) "calchim: Error; no O1D tracer !!!"
227	stop
228	else
229	nbq = nbq + 1
230	niq(nbq) = i_o1d
231	end if
232	i_o2 = igcm_o2
233	if (i_o2 == 0) then
234	write(,) "calchim: Error; no O2 tracer !!!"
235	stop
236	else
237	nbq = nbq + 1
238	niq(nbq) = i_o2
239	end if
240	i_o3 = igcm_o3
241	if (i_o3 == 0) then
242	write(,) "calchim: Error; no O3 tracer !!!"
243	stop
244	else
245	nbq = nbq + 1
246	niq(nbq) = i_o3
247	end if
248	i_h = igcm_h
249	if (i_h == 0) then
250	write(,) "calchim: Error; no H tracer !!!"
251	stop
252	else
253	nbq = nbq + 1
254	niq(nbq) = i_h
255	end if
256	i_h2 = igcm_h2
257	if (i_h2 == 0) then
258	write(,) "calchim: Error; no H2 tracer !!!"
259	stop
260	else
261	nbq = nbq + 1
262	niq(nbq) = i_h2
263	end if
264	i_oh = igcm_oh
265	if (i_oh == 0) then
266	write(,) "calchim: Error; no OH tracer !!!"
267	stop
268	else
269	nbq = nbq + 1
270	niq(nbq) = i_oh
271	end if
272	i_ho2 = igcm_ho2
273	if (i_ho2 == 0) then
274	write(,) "calchim: Error; no HO2 tracer !!!"
275	stop
276	else
277	nbq = nbq + 1
278	niq(nbq) = i_ho2
279	end if
280	i_h2o2 = igcm_h2o2
281	if (i_h2o2 == 0) then
282	write(,) "calchim: Error; no H2O2 tracer !!!"
283	stop
284	else
285	nbq = nbq + 1
286	niq(nbq) = i_h2o2
287	end if
288	i_ch4 = igcm_ch4
289	if (i_ch4 == 0) then
290	write(,) "calchim: Error; no CH4 tracer !!!"
291	write(,) "CH4 will be ignored in the chemistry"
292	else
293	nbq = nbq + 1
294	niq(nbq) = i_ch4
295	end if
296	i_n2 = igcm_n2
297	if (i_n2 == 0) then
298	write(,) "calchim: Error; no N2 tracer !!!"
299	stop
300	else
301	nbq = nbq + 1
302	niq(nbq) = i_n2
303	end if
304	i_n = igcm_n
305	if (i_n == 0) then
306	if (photochem) then
307	write(,) "calchim: Error; no N tracer !!!"
308	stop
309	end if
310	else
311	nbq = nbq + 1
312	niq(nbq) = i_n
313	end if
314	i_n2d = igcm_n2d
315	if (i_n2d == 0) then
316	if (photochem) then
317	write(,) "calchim: Error; no N2D tracer !!!"
318	stop
319	end if
320	else
321	nbq = nbq + 1
322	niq(nbq) = i_n2d
323	end if
324	i_no = igcm_no
325	if (i_no == 0) then
326	if (photochem) then
327	write(,) "calchim: Error; no NO tracer !!!"
328	stop
329	end if
330	else
331	nbq = nbq + 1
332	niq(nbq) = i_no
333	end if
334	i_no2 = igcm_no2
335	if (i_no2 == 0) then
336	if (photochem) then
337	write(,) "calchim: Error; no NO2 tracer !!!"
338	stop
339	end if
340	else
341	nbq = nbq + 1
342	niq(nbq) = i_no2
343	end if
344	i_h2o = igcm_h2o_vap
345	if (i_h2o == 0) then
346	write(,) "calchim: Error; no water vapor tracer !!!"
347	stop
348	else
349	nbq = nbq + 1
350	niq(nbq) = i_h2o
351	end if
352	!Check tracers needed for thermospheric chemistry
353	if(thermochem) then
354	chemthermod=0 !Default: C/O/H chemistry
355	!Nitrogen chemistry
356	!NO is used to determine if N chemistry is wanted
357	!chemthermod=2 -> N chemistry
358	if (i_no == 0) then
359	write(,) "calchim: no NO tracer"
360	write(,) "C/O/H themosp chemistry only "
361	else
362	chemthermod=2
363	write(,) "calchim: NO in traceur.def"
364	write(,) "Nitrogen chemistry included"
365	end if
366	! N
367	if(chemthermod == 2) then
368	if (i_n == 0) then
369	write(,) "calchim: Error; no N tracer !!!"
370	write(,) "N is needed if NO is in traceur.def"
371	stop
372	end if
373	! NO2
374	if (i_no2 == 0) then
375	write(,) "calchim: Error; no NO2 tracer !!!"
376	write(,) "NO2 is needed if NO is in traceur.def"
377	stop
378	end if
379	! N(2D)
380	if (i_n2d == 0) then
381	write(,) "calchim: Error; no N2D !!!"
382	write(,) "N2D is needed if NO is in traceur.def"
383	stop
384	end if
385	endif !Of if(chemthermod == 2)
386	! Ions
387	! O2+ is used to determine if ion chemistry is needed
388	! chemthermod=3 -> ion chemistry
389	i_o2plus = igcm_o2plus
390	if(chemthermod == 2) then
391	if (i_o2plus == 0) then
392	write(,) "calchim: no O2+ tracer; no ion chemistry"
393	else
394	nbq = nbq + 1
395	niq(nbq) = i_o2plus
396	chemthermod = 3
397	write(,) "calchim: O2+ in traceur.def"
398	write(,) "Ion chemistry included"
399	end if
400	else
401	if (i_o2plus /= 0) then
402	write(,) "calchim: O2+ is present, but NO is not!!!"
403	write(,) "Both must be in traceur.def if ionosphere wanted"
404	stop
405	endif
406	endif !Of if(chemthermod == 2)
407	! CO2+
408	i_co2plus = igcm_co2plus
409	if(chemthermod == 3) then
410	if (i_co2plus == 0) then
411	write(,) "calchim: Error; no CO2+ tracer !!!"
412	write(,) "CO2+ is needed if O2+ is in traceur.def"
413	stop
414	else
415	nbq = nbq + 1
416	niq(nbq) = i_co2plus
417	end if
418	else
419	if (i_co2plus /= 0) then
420	write(,) "calchim: Error: CO2+ is present, but O2+ is not!!!"
421	write(,) "Both must be in traceur.def if ionosphere wanted"
422	stop
423	endif
424	endif !Of if(chemthermod == 3)
425	! O+
426	i_oplus = igcm_oplus
427	if(chemthermod == 3) then
428	if (i_oplus == 0) then
429	write(,) "calchim: Error; no O+ tracer !!!"
430	write(,) "O+ is needed if O2+ is in traceur.def"
431	stop
432	else
433	nbq = nbq + 1
434	niq(nbq) = i_oplus
435	end if
436	else
437	if (i_oplus /= 0) then
438	write(,) "calchim: Error: O+ is present, but O2+ is not!!!"
439	write(,) "Both must be in traceur.def if ionosphere wanted"
440	stop
441	endif
442	endif !Of if (chemthermod == 3)
443	! CO+
444	i_coplus = igcm_coplus
445	if(chemthermod == 3) then
446	if (i_coplus == 0) then
447	write(,) "calchim: Error; no CO+ tracer !!!"
448	write(,) "CO+ is needed if O2+ is in traceur.def"
449	stop
450	else
451	nbq = nbq + 1
452	niq(nbq) = i_coplus
453	end if
454	else
455	if (i_coplus /= 0) then
456	write(,) "calchim: Error: CO+ is present, but O2+ is not!!!"
457	write(,) " Both must be in traceur.def if ionosphere wanted"
458	stop
459	endif
460	endif ! Of if (chemthermod == 3)
461	! C+
462	i_cplus = igcm_cplus
463	if(chemthermod == 3) then
464	if (i_cplus == 0) then
465	write(,) "calchim: Error; no C+ tracer !!!"
466	write(,) "C+ is needed if O2+ is in traceur.def"
467	stop
468	else
469	nbq = nbq + 1
470	niq(nbq) = i_cplus
471	end if
472	else
473	if (i_cplus /= 0) then
474	write(,) "calchim: Error; C+ is present, but O2+ is not!!!"
475	write(,) "Both must be in traceur.def if ionosphere wanted"
476	stop
477	endif
478	endif ! Of if (chemthermod == 3)
479	! N+
480	i_nplus = igcm_nplus
481	if(chemthermod == 3) then
482	if (i_nplus == 0) then
483	write(,) "calchim: Error; no N+ tracer !!!"
484	write(,) "N+ is needed if O2+ is in traceur.def"
485	stop
486	else
487	nbq = nbq + 1
488	niq(nbq) = i_nplus
489	end if
490	else
491	if (i_nplus /= 0) then
492	write(,) "calchim: Error: N+ is present, but O2+ is not!!!"
493	write(,) "Both must be in traceur.def if ionosphere wanted"
494	stop
495	endif
496	endif !Of if (chemthermod == 3)
497	! NO+
498	i_noplus = igcm_noplus
499	if(chemthermod == 3) then
500	if (i_noplus == 0) then
501	write(,) "calchim: Error; no NO+ tracer !!!"
502	write(,) "NO+ is needed if O2+ is in traceur.def"
503	stop
504	else
505	nbq = nbq + 1
506	niq(nbq) = i_noplus
507	end if
508	else
509	if (i_noplus /= 0) then
510	write(,) "calchim: Error: NO+ is present, but O2+ is not!!!"
511	write(,) "Both must be in traceur.def if ionosphere wanted"
512	stop
513	endif
514	endif !Of if (chemthermod == 3)
515	! N2+
516	i_n2plus = igcm_n2plus
517	if (chemthermod == 3) then
518	if (i_n2plus == 0) then
519	write(,) "calchim: Error; no N2+ tracer !!!"
520	write(,) "N2+ is needed if O2+ is in traceur.def"
521	stop
522	else
523	nbq = nbq + 1
524	niq(nbq) = i_n2plus
525	end if
526	else
527	if (i_n2plus /= 0) then
528	write(,) "calchim: Error: N2+ is present, but O2+ is not!!!"
529	write(,) "Both must be in traceur.def if ionosphere wanted"
530	stop
531	endif
532	endif !Of if (chemthermod == 3)
533	!H+
534	i_hplus = igcm_hplus
535	if (chemthermod == 3) then
536	if (i_hplus == 0) then
537	write(,) "calchim: Error; no H+ tracer !!!"
538	write(,) "H+ is needed if O2+ is in traceur.def"
539	stop
540	else
541	nbq = nbq + 1
542	niq(nbq) = i_hplus
543	end if
544	else
545	if (i_hplus /= 0) then
546	write(,) "calchim: Error: H+ is present, but O2+ is not!!!"
547	write(,) "Both must be in traceur.def if ionosphere wanted"
548	stop
549	endif
550	endif !Of if (chemthermod == 3)
551	! HCO2+
552	i_hco2plus = igcm_hco2plus
553	if(chemthermod == 3) then
554	if (i_hco2plus == 0) then
555	write(,) "calchim: Error; no HCO2+ tracer !!!"
556	write(,) "HCO2+ is needed if O2+ is in traceur.def"
557	stop
558	else
559	nbq = nbq + 1
560	niq(nbq) = i_hco2plus
561	end if
562	else
563	if (i_hco2plus /= 0) then
564	write(,) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
565	write(,) "Both must be in traceur.def if ionosphere wanted"
566	stop
567	endif
568	endif !Of if(chemthermod == 3)
569	!e-
570	i_elec = igcm_elec
571	if(chemthermod == 3) then
572	if (i_elec == 0) then
573	write(,) "calchim: Error; no e- tracer !!!"
574	write(,) "e- is needed if O2+ is in traceur.def"
575	stop
576	else
577	nbq = nbq + 1
578	niq(nbq) = i_elec
579	end if
580	else
581	if(i_elec /= 0) then
582	write(,) "calchim: Error: e- is present, but O2+ is not!!!"
583	write(,) "Both must be in traceur.def if ionosphere wanted"
584	stop
585	endif
586	endif !Of if (chemthermod == 3)
587	endif !Of thermochem
588
589	write(,) 'calchim: found nbq = ',nbq,' tracers'
590
591	firstcall = .false.
592	end if ! if (firstcall)
593
594	! Initializations
595
596	zycol(:,:) = 0.
597	dqchim(:,:,:) = 0.
598	dqschim(:,:) = 0.
599
600	! latvl1= 22.27
601	! lonvl1= -47.94
602	! ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim + &
603	! int(1.5+(lonvl1+180)*iim/360.)
604
605	!=======================================================================
606	! loop over grid
607	!=======================================================================
608
609	do ig = 1,ngrid
610
611	foundswitch = 0
612	do l = 1,nlayer
613	do i = 1,nbq
614	iq = niq(i) ! get tracer index
615	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
616	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
617	end do
618	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
619	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
620	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
621	zpress(l) = pplay(ig,l)/100.
622	ztemp(l) = zt(ig,l)
623	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
624	zlocal(l) = zzlay(ig,l)/1000.
625	zmmean(l) = mmean(ig,l)
626
627	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
628
629	surfdust1d(l) = surfdust(ig,l)*1.e-2
630	surfice1d(l) = surfice(ig,l)*1.e-2
631
632	! search for switch index between regions
633
634	if (photochem .and. thermochem) then
635	if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
636	lswitch = l
637	foundswitch = 1
638	end if
639	end if
640	if (.not. photochem) then
641	lswitch = 22
642	end if
643	if (.not. thermochem) then
644	lswitch = min(50,nlayer+1)
645	end if
646
647	end do ! of do l=1,nlayer
648
649	szacol = acos(mu0(ig))*180./pi
650	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
651
652	!=======================================================================
653	! call chemical subroutines
654	!=======================================================================
655
656	! chemistry in lower atmosphere
657
658	if (photochem) then
659
660	call photochemistry_asis(nlayer,nq, &
661	ig,lswitch,zycol,szacol,ptimestep, &
662	zpress,ztemp,zdens,zmmean,dist_sol, &
663	surfdust1d,surfice1d,jo3,taucol,iter)
664
665	! ozone photolysis, for output
666
667	do l = 1,nlayer
668	jo3_3d(ig,l) = jo3(l)
669	iter_3d(ig,l) = iter(l)
670	end do
671
672	! condensation of h2o2
673
674	call perosat(ngrid, nlayer, nq, &
675	ig,ptimestep,pplev,pplay, &
676	ztemp,zycol,dqcloud,dqscloud)
677	end if
678
679	! chemistry in upper atmosphere
680
681	if (thermochem) then
682	call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
683	zdens,zpress,zlocal,szacol,ptimestep,zday,&
684	em_no,em_o2)
685	do l=1,nlayer
686	emission_no(ig,l)=em_no(l)
687	emission_o2(ig,l)=em_o2(l)
688	enddo
689	end if
690
691	! dry deposition
692
693	if (depos) then
694	call deposition(ngrid, nlayer, nq, &
695	ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
696	zu, zv, zt, zycol, ptimestep, co2ice)
697	end if
698
699	!=======================================================================
700	! tendencies
701	!=======================================================================
702
703	! index of the most abundant species at each level
704
705	! major(:) = maxloc(zycol, dim = 2)
706
707	! tendency for the most abundant species = - sum of others
708
709	do l = 1,nlayer
710	iloc=maxloc(zycol(l,:))
711	iqmax=iloc(1)
712	do i = 1,nbq
713	iq = niq(i) ! get tracer index
714	if (iq /= iqmax) then
715	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
716	- zq(ig,l,iq))/ptimestep
717	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
718	- dqchim(ig,l,iq)
719	end if
720	end do
721	end do ! of do l = 1,nlayer
722
723	!=======================================================================
724	! end of loop over grid
725	!=======================================================================
726
727	end do ! of do ig=1,ngrid
728
729	!=======================================================================
730	! write outputs
731	!=======================================================================
732
733	! value of parameter 'output' to trigger writting of outputs
734	! is set above at the declaration of the variable.
735
736	if (photochem .and. output) then
737	if (ngrid > 1) then
738	call writediagfi(ngrid,'jo3','j o3->o1d', &
739	's-1',3,jo3_3d(1,1))
740	call writediagfi(ngrid,'iter','iterations', &
741	' ',3,iter_3d(1,1))
742	call writediagfi(ngrid,'emission_no', &
743	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
744	call writediagfi(ngrid,'emission_o2', &
745	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
746
747	if (callstats) then
748	call wstats(ngrid,'jo3','j o3->o1d', &
749	's-1',3,jo3_3d(1,1))
750	call wstats(ngrid,'emission_no', &
751	'NO nightglow emission rate','cm-3 s-1',3,emission_no)
752	call wstats(ngrid,'emission_o2', &
753	'O2 nightglow emission rate','cm-3 s-1',3,emission_o2)
754	call wstats(ngrid,'mmean','mean molecular mass', &
755	'g.mole-1',3,mmean(1,1))
756	endif
757	end if ! of if (ngrid.gt.1)
758	end if ! of if (output)
759
760	return
761	end

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

Download in other formats:

Original Format