Context Navigation

calchim.F90 @ 2029

Last change on this file since 2029 was 2029, checked in by flefevre, 7 years ago
Suite du chantier photolyses on-line. Routines non operationnelles par defaut, controle des resultats en cours...
File size: 28.0 KB

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

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

Download in other formats:

Original Format