Context Navigation

calchim.F90 @ 658

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

Mars GCM:

updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h, paramfoto_compact.F, param_read.F , thermosphere.F).
minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function which seems to be a problem for g95) .
minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if none is computed.
in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP" and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines defined in "moldiff_red.F". Added use of automatic-sized array in "tridag" and "tridag_sp" local array "gam" (to avoid using an a priori oversized local array).

FGG+JYC+EM

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

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

Download in other formats:

Original Format