Context Navigation

calchim.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: 26.3 KB

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

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

Download in other formats:

Original Format