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