Context Navigation

calchim.F90 @ 1775

Last change on this file since 1775 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format