Context Navigation

calchim_asis.F90 @ 1541

Last change on this file since 1541 was 1495, checked in by flefevre, 10 years ago
Implementation du coeur chimique Adaptative Semi-Implicit Symmetric (ASIS) Par defaut non actif pour l'instant (variable logique asis = .false.)
File size: 26.9 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
170	integer :: iter(nlayer) ! Number of chemical iterations
171	! within one physical timestep
172
173	! for output:
174
175	logical :: output ! to issue calls to writediagfi and stats
176	parameter (output = .true.)
177	real :: jo3_3d(ngrid,nlayer) ! Photodissociation rate O3->O1D (s-1)
178	real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
179	! within one physical timestep
180
181	!=======================================================================
182	! initialization of the chemistry (first call only)
183	!=======================================================================
184
185	if (firstcall) then
186
187	if (photochem) then
188	print*,'calchim: Read photolysis lookup table'
189	call read_phototable
190	end if
191	! find index of chemical tracers to use
192	allocate(niq(nq))
193	! Listed here are all tracers that can go into photochemistry
194	nbq = 0 ! to count number of tracers
195	! Species ALWAYS present if photochem=.T. or thermochem=.T.
196	i_co2 = igcm_co2
197	if (i_co2 == 0) then
198	write(,) "calchim: Error; no CO2 tracer !!!"
199	stop
200	else
201	nbq = nbq + 1
202	niq(nbq) = i_co2
203	end if
204	i_co = igcm_co
205	if (i_co == 0) then
206	write(,) "calchim: Error; no CO tracer !!!"
207	stop
208	else
209	nbq = nbq + 1
210	niq(nbq) = i_co
211	end if
212	i_o = igcm_o
213	if (i_o == 0) then
214	write(,) "calchim: Error; no O tracer !!!"
215	stop
216	else
217	nbq = nbq + 1
218	niq(nbq) = i_o
219	end if
220	i_o1d = igcm_o1d
221	if (i_o1d == 0) then
222	write(,) "calchim: Error; no O1D tracer !!!"
223	stop
224	else
225	nbq = nbq + 1
226	niq(nbq) = i_o1d
227	end if
228	i_o2 = igcm_o2
229	if (i_o2 == 0) then
230	write(,) "calchim: Error; no O2 tracer !!!"
231	stop
232	else
233	nbq = nbq + 1
234	niq(nbq) = i_o2
235	end if
236	i_o3 = igcm_o3
237	if (i_o3 == 0) then
238	write(,) "calchim: Error; no O3 tracer !!!"
239	stop
240	else
241	nbq = nbq + 1
242	niq(nbq) = i_o3
243	end if
244	i_h = igcm_h
245	if (i_h == 0) then
246	write(,) "calchim: Error; no H tracer !!!"
247	stop
248	else
249	nbq = nbq + 1
250	niq(nbq) = i_h
251	end if
252	i_h2 = igcm_h2
253	if (i_h2 == 0) then
254	write(,) "calchim: Error; no H2 tracer !!!"
255	stop
256	else
257	nbq = nbq + 1
258	niq(nbq) = i_h2
259	end if
260	i_oh = igcm_oh
261	if (i_oh == 0) then
262	write(,) "calchim: Error; no OH tracer !!!"
263	stop
264	else
265	nbq = nbq + 1
266	niq(nbq) = i_oh
267	end if
268	i_ho2 = igcm_ho2
269	if (i_ho2 == 0) then
270	write(,) "calchim: Error; no HO2 tracer !!!"
271	stop
272	else
273	nbq = nbq + 1
274	niq(nbq) = i_ho2
275	end if
276	i_h2o2 = igcm_h2o2
277	if (i_h2o2 == 0) then
278	write(,) "calchim: Error; no H2O2 tracer !!!"
279	stop
280	else
281	nbq = nbq + 1
282	niq(nbq) = i_h2o2
283	end if
284	i_ch4 = igcm_ch4
285	if (i_ch4 == 0) then
286	write(,) "calchim: Error; no CH4 tracer !!!"
287	write(,) "CH4 will be ignored in the chemistry"
288	else
289	nbq = nbq + 1
290	niq(nbq) = i_ch4
291	end if
292	i_n2 = igcm_n2
293	if (i_n2 == 0) then
294	write(,) "calchim: Error; no N2 tracer !!!"
295	stop
296	else
297	nbq = nbq + 1
298	niq(nbq) = i_n2
299	end if
300	i_n = igcm_n
301	if (i_n == 0) then
302	if (photochem) then
303	write(,) "calchim: Error; no N tracer !!!"
304	stop
305	end if
306	else
307	nbq = nbq + 1
308	niq(nbq) = i_n
309	end if
310	i_n2d = igcm_n2d
311	if (i_n2d == 0) then
312	if (photochem) then
313	write(,) "calchim: Error; no N2D tracer !!!"
314	stop
315	end if
316	else
317	nbq = nbq + 1
318	niq(nbq) = i_n2d
319	end if
320	i_no = igcm_no
321	if (i_no == 0) then
322	if (photochem) then
323	write(,) "calchim: Error; no NO tracer !!!"
324	stop
325	end if
326	else
327	nbq = nbq + 1
328	niq(nbq) = i_no
329	end if
330	i_no2 = igcm_no2
331	if (i_no2 == 0) then
332	if (photochem) then
333	write(,) "calchim: Error; no NO2 tracer !!!"
334	stop
335	end if
336	else
337	nbq = nbq + 1
338	niq(nbq) = i_no2
339	end if
340	i_h2o = igcm_h2o_vap
341	if (i_h2o == 0) then
342	write(,) "calchim: Error; no water vapor tracer !!!"
343	stop
344	else
345	nbq = nbq + 1
346	niq(nbq) = i_h2o
347	end if
348	!Check tracers needed for thermospheric chemistry
349	if(thermochem) then
350	chemthermod=0 !Default: C/O/H chemistry
351	!Nitrogen chemistry
352	!NO is used to determine if N chemistry is wanted
353	!chemthermod=2 -> N chemistry
354	if (i_no == 0) then
355	write(,) "calchim: no NO tracer"
356	write(,) "C/O/H themosp chemistry only "
357	else
358	chemthermod=2
359	write(,) "calchim: NO in traceur.def"
360	write(,) "Nitrogen chemistry included"
361	end if
362	! N
363	if(chemthermod == 2) then
364	if (i_n == 0) then
365	write(,) "calchim: Error; no N tracer !!!"
366	write(,) "N is needed if NO is in traceur.def"
367	stop
368	end if
369	! NO2
370	if (i_no2 == 0) then
371	write(,) "calchim: Error; no NO2 tracer !!!"
372	write(,) "NO2 is needed if NO is in traceur.def"
373	stop
374	end if
375	! N(2D)
376	if (i_n2d == 0) then
377	write(,) "calchim: Error; no N2D !!!"
378	write(,) "N2D is needed if NO is in traceur.def"
379	stop
380	end if
381	endif !Of if(chemthermod == 2)
382	! Ions
383	! O2+ is used to determine if ion chemistry is needed
384	! chemthermod=3 -> ion chemistry
385	i_o2plus = igcm_o2plus
386	if(chemthermod == 2) then
387	if (i_o2plus == 0) then
388	write(,) "calchim: no O2+ tracer; no ion chemistry"
389	else
390	nbq = nbq + 1
391	niq(nbq) = i_o2plus
392	chemthermod = 3
393	write(,) "calchim: O2+ in traceur.def"
394	write(,) "Ion chemistry included"
395	end if
396	else
397	if (i_o2plus /= 0) then
398	write(,) "calchim: O2+ is present, but NO is not!!!"
399	write(,) "Both must be in traceur.def if ionosphere wanted"
400	stop
401	endif
402	endif !Of if(chemthermod == 2)
403	! CO2+
404	i_co2plus = igcm_co2plus
405	if(chemthermod == 3) then
406	if (i_co2plus == 0) then
407	write(,) "calchim: Error; no CO2+ tracer !!!"
408	write(,) "CO2+ is needed if O2+ is in traceur.def"
409	stop
410	else
411	nbq = nbq + 1
412	niq(nbq) = i_co2plus
413	end if
414	else
415	if (i_co2plus /= 0) then
416	write(,) "calchim: Error: CO2+ is present, but O2+ is not!!!"
417	write(,) "Both must be in traceur.def if ionosphere wanted"
418	stop
419	endif
420	endif !Of if(chemthermod == 3)
421	! O+
422	i_oplus = igcm_oplus
423	if(chemthermod == 3) then
424	if (i_oplus == 0) then
425	write(,) "calchim: Error; no O+ tracer !!!"
426	write(,) "O+ is needed if O2+ is in traceur.def"
427	stop
428	else
429	nbq = nbq + 1
430	niq(nbq) = i_oplus
431	end if
432	else
433	if (i_oplus /= 0) then
434	write(,) "calchim: Error: O+ is present, but O2+ is not!!!"
435	write(,) "Both must be in traceur.def if ionosphere wanted"
436	stop
437	endif
438	endif !Of if (chemthermod == 3)
439	! CO+
440	i_coplus = igcm_coplus
441	if(chemthermod == 3) then
442	if (i_coplus == 0) then
443	write(,) "calchim: Error; no CO+ tracer !!!"
444	write(,) "CO+ is needed if O2+ is in traceur.def"
445	stop
446	else
447	nbq = nbq + 1
448	niq(nbq) = i_coplus
449	end if
450	else
451	if (i_coplus /= 0) then
452	write(,) "calchim: Error: CO+ is present, but O2+ is not!!!"
453	write(,) " Both must be in traceur.def if ionosphere wanted"
454	stop
455	endif
456	endif ! Of if (chemthermod == 3)
457	! C+
458	i_cplus = igcm_cplus
459	if(chemthermod == 3) then
460	if (i_cplus == 0) then
461	write(,) "calchim: Error; no C+ tracer !!!"
462	write(,) "C+ is needed if O2+ is in traceur.def"
463	stop
464	else
465	nbq = nbq + 1
466	niq(nbq) = i_cplus
467	end if
468	else
469	if (i_cplus /= 0) then
470	write(,) "calchim: Error; C+ is present, but O2+ is not!!!"
471	write(,) "Both must be in traceur.def if ionosphere wanted"
472	stop
473	endif
474	endif ! Of if (chemthermod == 3)
475	! N+
476	i_nplus = igcm_nplus
477	if(chemthermod == 3) then
478	if (i_nplus == 0) then
479	write(,) "calchim: Error; no N+ tracer !!!"
480	write(,) "N+ is needed if O2+ is in traceur.def"
481	stop
482	else
483	nbq = nbq + 1
484	niq(nbq) = i_nplus
485	end if
486	else
487	if (i_nplus /= 0) then
488	write(,) "calchim: Error: N+ is present, but O2+ is not!!!"
489	write(,) "Both must be in traceur.def if ionosphere wanted"
490	stop
491	endif
492	endif !Of if (chemthermod == 3)
493	! NO+
494	i_noplus = igcm_noplus
495	if(chemthermod == 3) then
496	if (i_noplus == 0) then
497	write(,) "calchim: Error; no NO+ tracer !!!"
498	write(,) "NO+ is needed if O2+ is in traceur.def"
499	stop
500	else
501	nbq = nbq + 1
502	niq(nbq) = i_noplus
503	end if
504	else
505	if (i_noplus /= 0) then
506	write(,) "calchim: Error: NO+ is present, but O2+ is not!!!"
507	write(,) "Both must be in traceur.def if ionosphere wanted"
508	stop
509	endif
510	endif !Of if (chemthermod == 3)
511	! N2+
512	i_n2plus = igcm_n2plus
513	if (chemthermod == 3) then
514	if (i_n2plus == 0) then
515	write(,) "calchim: Error; no N2+ tracer !!!"
516	write(,) "N2+ is needed if O2+ is in traceur.def"
517	stop
518	else
519	nbq = nbq + 1
520	niq(nbq) = i_n2plus
521	end if
522	else
523	if (i_n2plus /= 0) then
524	write(,) "calchim: Error: N2+ is present, but O2+ is not!!!"
525	write(,) "Both must be in traceur.def if ionosphere wanted"
526	stop
527	endif
528	endif !Of if (chemthermod == 3)
529	!H+
530	i_hplus = igcm_hplus
531	if (chemthermod == 3) then
532	if (i_hplus == 0) then
533	write(,) "calchim: Error; no H+ tracer !!!"
534	write(,) "H+ is needed if O2+ is in traceur.def"
535	stop
536	else
537	nbq = nbq + 1
538	niq(nbq) = i_hplus
539	end if
540	else
541	if (i_hplus /= 0) then
542	write(,) "calchim: Error: H+ is present, but O2+ is not!!!"
543	write(,) "Both must be in traceur.def if ionosphere wanted"
544	stop
545	endif
546	endif !Of if (chemthermod == 3)
547	! HCO2+
548	i_hco2plus = igcm_hco2plus
549	if(chemthermod == 3) then
550	if (i_hco2plus == 0) then
551	write(,) "calchim: Error; no HCO2+ tracer !!!"
552	write(,) "HCO2+ is needed if O2+ is in traceur.def"
553	stop
554	else
555	nbq = nbq + 1
556	niq(nbq) = i_hco2plus
557	end if
558	else
559	if (i_hco2plus /= 0) then
560	write(,) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
561	write(,) "Both must be in traceur.def if ionosphere wanted"
562	stop
563	endif
564	endif !Of if(chemthermod == 3)
565	!e-
566	i_elec = igcm_elec
567	if(chemthermod == 3) then
568	if (i_elec == 0) then
569	write(,) "calchim: Error; no e- tracer !!!"
570	write(,) "e- is needed if O2+ is in traceur.def"
571	stop
572	else
573	nbq = nbq + 1
574	niq(nbq) = i_elec
575	end if
576	else
577	if(i_elec /= 0) then
578	write(,) "calchim: Error: e- is present, but O2+ is not!!!"
579	write(,) "Both must be in traceur.def if ionosphere wanted"
580	stop
581	endif
582	endif !Of if (chemthermod == 3)
583	endif !Of thermochem
584
585	write(,) 'calchim: found nbq = ',nbq,' tracers'
586
587	firstcall = .false.
588	end if ! if (firstcall)
589
590	! Initializations
591
592	zycol(:,:) = 0.
593	dqchim(:,:,:) = 0.
594	dqschim(:,:) = 0.
595
596	! latvl1= 22.27
597	! lonvl1= -47.94
598	! ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim + &
599	! int(1.5+(lonvl1+180)*iim/360.)
600
601	!=======================================================================
602	! loop over grid
603	!=======================================================================
604
605	do ig = 1,ngrid
606
607	foundswitch = 0
608	do l = 1,nlayer
609	do i = 1,nbq
610	iq = niq(i) ! get tracer index
611	zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
612	zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
613	end do
614	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
615	zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
616	zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
617	zpress(l) = pplay(ig,l)/100.
618	ztemp(l) = zt(ig,l)
619	zdens(l) = zpress(l)/(kb1.e4ztemp(l))
620	zlocal(l) = zzlay(ig,l)/1000.
621	zmmean(l) = mmean(ig,l)
622
623	! surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
624
625	surfdust1d(l) = surfdust(ig,l)*1.e-2
626	surfice1d(l) = surfice(ig,l)*1.e-2
627
628	! search for switch index between regions
629
630	if (photochem .and. thermochem) then
631	if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
632	lswitch = l
633	foundswitch = 1
634	end if
635	end if
636	if (.not. photochem) then
637	lswitch = 22
638	end if
639	if (.not. thermochem) then
640	lswitch = min(50,nlayer+1)
641	end if
642
643	end do ! of do l=1,nlayer
644
645	szacol = acos(mu0(ig))*180./pi
646	taucol = tauref(ig)*(700./610.) ! provisoire en attente de nouveau jmars
647
648	!=======================================================================
649	! call chemical subroutines
650	!=======================================================================
651
652	! chemistry in lower atmosphere
653
654	if (photochem) then
655
656	call photochemistry_asis(nlayer,nq, &
657	ig,lswitch,zycol,szacol,ptimestep, &
658	zpress,ztemp,zdens,zmmean,dist_sol, &
659	surfdust1d,surfice1d,jo3,taucol,iter)
660
661	! ozone photolysis, for output
662
663	do l = 1,nlayer
664	jo3_3d(ig,l) = jo3(l)
665	iter_3d(ig,l) = iter(l)
666	end do
667
668	! condensation of h2o2
669
670	call perosat(ngrid, nlayer, nq, &
671	ig,ptimestep,pplev,pplay, &
672	ztemp,zycol,dqcloud,dqscloud)
673	end if
674
675	! chemistry in upper atmosphere
676
677	if (thermochem) then
678	call chemthermos(ig,nlayer,lswitch,chemthermod,zycol,ztemp,&
679	zdens,zpress,zlocal,szacol,ptimestep,zday)
680	end if
681
682	! dry deposition
683
684	if (depos) then
685	call deposition(ngrid, nlayer, nq, &
686	ig, ig_vl1, pplay, pplev, zzlay, zzlev, &
687	zu, zv, zt, zycol, ptimestep, co2ice)
688	end if
689
690	!=======================================================================
691	! tendencies
692	!=======================================================================
693
694	! index of the most abundant species at each level
695
696	! major(:) = maxloc(zycol, dim = 2)
697
698	! tendency for the most abundant species = - sum of others
699
700	do l = 1,nlayer
701	iloc=maxloc(zycol(l,:))
702	iqmax=iloc(1)
703	do i = 1,nbq
704	iq = niq(i) ! get tracer index
705	if (iq /= iqmax) then
706	dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) &
707	- zq(ig,l,iq))/ptimestep
708	dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) &
709	- dqchim(ig,l,iq)
710	end if
711	end do
712	end do ! of do l = 1,nlayer
713
714	!=======================================================================
715	! end of loop over grid
716	!=======================================================================
717
718	end do ! of do ig=1,ngrid
719
720	!=======================================================================
721	! write outputs
722	!=======================================================================
723
724	! value of parameter 'output' to trigger writting of outputs
725	! is set above at the declaration of the variable.
726
727	if (photochem .and. output) then
728	if (ngrid > 1) then
729	call writediagfi(ngrid,'jo3','j o3->o1d', &
730	's-1',3,jo3_3d(1,1))
731	call writediagfi(ngrid,'iter','iterations', &
732	' ',3,iter_3d(1,1))
733	if (callstats) then
734	call wstats(ngrid,'jo3','j o3->o1d', &
735	's-1',3,jo3_3d(1,1))
736	call wstats(ngrid,'mmean','mean molecular mass', &
737	'g.mole-1',3,mmean(1,1))
738	endif
739	end if ! of if (ngrid.gt.1)
740	end if ! of if (output)
741
742	return
743	end

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

Download in other formats:

Original Format