Context Navigation

calchim.F90 @ 635

Last change on this file since 635 was 635, checked in by emillour, 13 years ago

Mars GCM: Update of the chemistry package, including:

93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
computation of photoabsorption using raytracing
improved time stepping in the photochemistry
updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
Upadted surfacearea.F
Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM

File size: 25.7 KB

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

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

Download in other formats:

Original Format