1 | !========================================================================== |
---|
2 | |
---|
3 | subroutine photolysis_asis(nlayer, & |
---|
4 | lswitch, press, temp, sza, tauref, & |
---|
5 | zmmean, dist_sol, rmco2, rmo3, v_phot) |
---|
6 | |
---|
7 | !========================================================================== |
---|
8 | |
---|
9 | use comcstfi_h |
---|
10 | |
---|
11 | implicit none |
---|
12 | |
---|
13 | #include "chimiedata.h" |
---|
14 | |
---|
15 | !========================================================================== |
---|
16 | ! input: |
---|
17 | !========================================================================== |
---|
18 | |
---|
19 | integer, intent(in) :: nlayer ! number of atmospheric layers |
---|
20 | integer :: lswitch ! interface level between chemistries |
---|
21 | real :: press(nlayer) ! pressure (hPa) |
---|
22 | real :: temp(nlayer) ! temperature (K) |
---|
23 | real :: sza ! solar zenith angle (deg) |
---|
24 | real :: tauref ! optical depth at 7 hpa |
---|
25 | real :: zmmean(nlayer) ! mean molecular mass (g) |
---|
26 | real :: dist_sol ! sun distance (AU) |
---|
27 | real :: rmco2(nlayer) ! co2 mixing ratio |
---|
28 | real :: rmo3(nlayer) ! ozone mixing ratio |
---|
29 | |
---|
30 | !========================================================================== |
---|
31 | ! output: interpolated photodissociation rates (s-1) |
---|
32 | !========================================================================== |
---|
33 | |
---|
34 | real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot |
---|
35 | |
---|
36 | !========================================================================== |
---|
37 | ! local: |
---|
38 | !========================================================================== |
---|
39 | |
---|
40 | integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp, & |
---|
41 | iozo, isza, itau, it, l |
---|
42 | |
---|
43 | integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, & |
---|
44 | j_o3_o, j_h2o, j_hdo, j_h2o2, j_ho2, j_no, j_no2, & |
---|
45 | j_hno3, j_hno4, & |
---|
46 | j_ch4_ch3_h, j_ch4_1ch2_h2, j_ch4_3ch2_h_h, & |
---|
47 | j_ch4_ch_h2_h, j_ch3o2h, j_ch2o_hco, j_ch2o_co, & |
---|
48 | j_ch3oh, j_c2h6, j_hcl, j_hocl, j_clo, j_so2, j_so, & |
---|
49 | j_h2s, j_so3 |
---|
50 | |
---|
51 | real :: col(nlayer) ! overhead air column (molecule cm-2) |
---|
52 | real :: colo3(nlayer) ! overhead ozone column (molecule cm-2) |
---|
53 | real :: poids(2,2,2,2,2) ! 5D interpolation weights |
---|
54 | real :: tref ! temperature at 1.9 hPa in the gcm (K) |
---|
55 | real :: table_temp(ntemp) ! temperatures at 1.9 hPa in jmars (K) |
---|
56 | real :: cinf, csup, cicol, ciozo, cisza, citemp, citau |
---|
57 | real :: colo3min, dp, coef |
---|
58 | real :: ratio_o3(nlayer) |
---|
59 | real :: tau |
---|
60 | real :: j(nlayer,nd) |
---|
61 | |
---|
62 | !========================================================================== |
---|
63 | ! day/night criterion |
---|
64 | !========================================================================== |
---|
65 | |
---|
66 | if (sza <= 95.) then |
---|
67 | |
---|
68 | !========================================================================== |
---|
69 | ! temperatures at 1.9 hPa in lookup table |
---|
70 | !========================================================================== |
---|
71 | |
---|
72 | table_temp(1) = 226.2 |
---|
73 | table_temp(2) = 206.2 |
---|
74 | table_temp(3) = 186.2 |
---|
75 | table_temp(4) = 169.8 |
---|
76 | |
---|
77 | !========================================================================== |
---|
78 | ! interpolation in solar zenith angle |
---|
79 | !========================================================================== |
---|
80 | |
---|
81 | indsza = nsza - 1 |
---|
82 | do isza = 1,nsza |
---|
83 | if (szatab(isza) >= sza) then |
---|
84 | indsza = isza - 1 |
---|
85 | indsza = max(indsza, 1) |
---|
86 | exit |
---|
87 | end if |
---|
88 | end do |
---|
89 | cisza = (sza - szatab(indsza)) & |
---|
90 | /(szatab(indsza + 1) - szatab(indsza)) |
---|
91 | |
---|
92 | !========================================================================== |
---|
93 | ! interpolation in dust (tau) |
---|
94 | !========================================================================== |
---|
95 | |
---|
96 | tau = min(tauref, tautab(ntau)) |
---|
97 | tau = max(tau, tautab(1)) |
---|
98 | |
---|
99 | indtau = ntau - 1 |
---|
100 | do itau = 1,ntau |
---|
101 | if (tautab(itau) >= tau) then |
---|
102 | indtau = itau - 1 |
---|
103 | indtau = max(indtau, 1) |
---|
104 | exit |
---|
105 | end if |
---|
106 | end do |
---|
107 | citau = (tau - tautab(indtau)) & |
---|
108 | /(tautab(indtau + 1) - tautab(indtau)) |
---|
109 | |
---|
110 | !========================================================================== |
---|
111 | ! co2 and ozone columns |
---|
112 | !========================================================================== |
---|
113 | |
---|
114 | ! co2 column at model top (molecule.cm-2) |
---|
115 | |
---|
116 | col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100. & |
---|
117 | /(zmmean(lswitch-1)*g) |
---|
118 | |
---|
119 | ! ozone column at model top |
---|
120 | |
---|
121 | colo3(lswitch-1) = 0. |
---|
122 | |
---|
123 | ! co2 and ozone columns for other levels (molecule.cm-2) |
---|
124 | |
---|
125 | do l = lswitch-2,1,-1 |
---|
126 | dp = (press(l) - press(l+1))*100. |
---|
127 | col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5 & |
---|
128 | *6.022e22*dp/(zmmean(l)*g) |
---|
129 | col(l) = min(col(l), colairtab(1)) |
---|
130 | colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 & |
---|
131 | *6.022e22*dp/(zmmean(l)*g) |
---|
132 | end do |
---|
133 | |
---|
134 | ! ratio ozone column/minimal theoretical column (0.1 micron-atm) |
---|
135 | |
---|
136 | ! ro3 = 7.171e-10 is the o3 mixing ratio for a uniform |
---|
137 | ! profile giving a column 0.1 micron-atmosphere at |
---|
138 | ! a surface pressure of 10 hpa. |
---|
139 | |
---|
140 | do l = 1,lswitch-1 |
---|
141 | colo3min = col(l)*7.171e-10 |
---|
142 | ratio_o3(l) = colo3(l)/colo3min |
---|
143 | ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.) |
---|
144 | ratio_o3(l) = max(ratio_o3(l), 1.) |
---|
145 | end do |
---|
146 | |
---|
147 | !========================================================================== |
---|
148 | ! temperature dependence |
---|
149 | !========================================================================== |
---|
150 | |
---|
151 | ! 1) search for temperature at 1.9 hPa (tref): vertical interpolation |
---|
152 | |
---|
153 | tref = temp(1) |
---|
154 | do l = (lswitch-1)-1,1,-1 |
---|
155 | if (press(l) > 1.9) then |
---|
156 | cinf = (press(l) - 1.9)/(press(l) - press(l+1)) |
---|
157 | csup = 1. - cinf |
---|
158 | tref = cinf*temp(l+1) + csup*temp(l) |
---|
159 | exit |
---|
160 | end if |
---|
161 | end do |
---|
162 | |
---|
163 | ! 2) interpolation in temperature |
---|
164 | |
---|
165 | tref = min(tref,table_temp(1)) |
---|
166 | tref = max(tref,table_temp(ntemp)) |
---|
167 | |
---|
168 | do it = 2, ntemp |
---|
169 | if (table_temp(it) <= tref) then |
---|
170 | citemp = (log(tref) - log(table_temp(it))) & |
---|
171 | /(log(table_temp(it-1)) - log(table_temp(it))) |
---|
172 | indtemp = it - 1 |
---|
173 | exit |
---|
174 | end if |
---|
175 | end do |
---|
176 | |
---|
177 | !========================================================================== |
---|
178 | ! loop over vertical levels |
---|
179 | !========================================================================== |
---|
180 | |
---|
181 | do l = 1,lswitch-1 |
---|
182 | |
---|
183 | ! interpolation in air column |
---|
184 | |
---|
185 | indcol = nz - 1 |
---|
186 | do icol = 1,nz |
---|
187 | if (colairtab(icol) < col(l)) then |
---|
188 | indcol = icol - 1 |
---|
189 | exit |
---|
190 | end if |
---|
191 | end do |
---|
192 | cicol = (log(col(l)) - log(colairtab(indcol + 1))) & |
---|
193 | /(log(colairtab(indcol)) - log(colairtab(indcol + 1))) |
---|
194 | |
---|
195 | ! interpolation in ozone column |
---|
196 | |
---|
197 | indozo = nozo - 1 |
---|
198 | do iozo = 1,nozo |
---|
199 | if (table_ozo(iozo)*10. >= ratio_o3(l)) then |
---|
200 | indozo = iozo - 1 |
---|
201 | indozo = max(indozo, 1) |
---|
202 | exit |
---|
203 | end if |
---|
204 | end do |
---|
205 | ciozo = (ratio_o3(l) - table_ozo(indozo)*10.) & |
---|
206 | /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.) |
---|
207 | |
---|
208 | ! 4-dimensional interpolation weights |
---|
209 | |
---|
210 | ! poids(temp,sza,co2,o3,tau) |
---|
211 | |
---|
212 | poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau) |
---|
213 | poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau) |
---|
214 | poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau) |
---|
215 | poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau) |
---|
216 | poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau) |
---|
217 | poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau) |
---|
218 | poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau) |
---|
219 | poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau) |
---|
220 | poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau) |
---|
221 | poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau) |
---|
222 | poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau) |
---|
223 | poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau) |
---|
224 | poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau) |
---|
225 | poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau) |
---|
226 | poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau) |
---|
227 | poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau) |
---|
228 | ! |
---|
229 | poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau |
---|
230 | poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau |
---|
231 | poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau |
---|
232 | poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau |
---|
233 | poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau |
---|
234 | poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau |
---|
235 | poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau |
---|
236 | poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau |
---|
237 | poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau |
---|
238 | poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau |
---|
239 | poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau |
---|
240 | poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau |
---|
241 | poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau |
---|
242 | poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau |
---|
243 | poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau |
---|
244 | poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau |
---|
245 | |
---|
246 | ! 4-dimensional interpolation in the lookup table |
---|
247 | |
---|
248 | do ij = 1,nd |
---|
249 | j(l,ij) = & |
---|
250 | poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij) & |
---|
251 | + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij) & |
---|
252 | + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij) & |
---|
253 | + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij) & |
---|
254 | + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij) & |
---|
255 | + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij) & |
---|
256 | + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij) & |
---|
257 | + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij) & |
---|
258 | + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij) & |
---|
259 | + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij) & |
---|
260 | + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij) & |
---|
261 | + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij) & |
---|
262 | + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij) & |
---|
263 | + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij) & |
---|
264 | + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij) & |
---|
265 | + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij) & |
---|
266 | ! |
---|
267 | + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij) & |
---|
268 | + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij) & |
---|
269 | + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij) & |
---|
270 | + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij) & |
---|
271 | + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij) & |
---|
272 | + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij) & |
---|
273 | + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij) & |
---|
274 | + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij) & |
---|
275 | + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij) & |
---|
276 | + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij) & |
---|
277 | + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij) & |
---|
278 | + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij) & |
---|
279 | + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij) & |
---|
280 | + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij) & |
---|
281 | + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij) & |
---|
282 | + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij) |
---|
283 | end do |
---|
284 | |
---|
285 | ! correction for sun distance |
---|
286 | |
---|
287 | do ij = 1,nd |
---|
288 | j(l,ij) = j(l,ij)*(1.52/dist_sol)**2. |
---|
289 | end do |
---|
290 | |
---|
291 | !========================================================================== |
---|
292 | ! end of loop over vertical levels |
---|
293 | !========================================================================== |
---|
294 | |
---|
295 | end do |
---|
296 | |
---|
297 | else |
---|
298 | |
---|
299 | !========================================================================== |
---|
300 | ! night |
---|
301 | !========================================================================== |
---|
302 | |
---|
303 | j(:,:) = 0. |
---|
304 | |
---|
305 | end if |
---|
306 | |
---|
307 | ! photodissociation rates numbering in the lookup table |
---|
308 | |
---|
309 | ! jmars.20140930 |
---|
310 | |
---|
311 | j_o2_o = 1 ! o2 + hv -> o + o |
---|
312 | j_o2_o1d = 2 ! o2 + hv -> o + o(1d) |
---|
313 | j_co2_o = 3 ! co2 + hv -> co + o |
---|
314 | j_co2_o1d = 4 ! co2 + hv -> co + o(1d) |
---|
315 | j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d) |
---|
316 | j_o3_o = 6 ! o3 + hv -> o2 + o |
---|
317 | j_h2o = 7 ! h2o + hv -> h + oh |
---|
318 | j_h2o2 = 8 ! h2o2 + hv -> oh + oh |
---|
319 | j_ho2 = 9 ! ho2 + hv -> oh + o |
---|
320 | j_no = 10 ! no + hv -> n + o |
---|
321 | j_no2 = 11 ! no2 + hv -> no + o |
---|
322 | j_hno3 = 12 ! hno3 + hv -> no2 + oh |
---|
323 | j_hno4 = 13 ! hno4 + hv -> no2 + ho2 |
---|
324 | |
---|
325 | ! jmars.20111014 |
---|
326 | |
---|
327 | ! j_o2_o = 1 ! o2 + hv -> o + o |
---|
328 | ! j_o2_o1d = 2 ! o2 + hv -> o + o(1d) |
---|
329 | ! j_co2_o = 3 ! co2 + hv -> co + o |
---|
330 | ! j_co2_o1d = 4 ! co2 + hv -> co + o(1d) |
---|
331 | ! j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d) |
---|
332 | ! j_o3_o = 6 ! o3 + hv -> o2 + o |
---|
333 | ! j_h2o = 7 ! h2o + hv -> h + oh |
---|
334 | ! j_hdo = 8 ! hdo + hv -> d + oh |
---|
335 | ! j_h2o2 = 9 ! h2o2 + hv -> oh + oh |
---|
336 | ! j_ho2 = 10 ! ho2 + hv -> oh + o |
---|
337 | ! j_no2 = 11 ! no2 + hv -> no + o |
---|
338 | ! j_ch4_ch3_h = 12 ! ch4 + hv -> ch3 + h |
---|
339 | ! j_ch4_1ch2_h2 = 13 ! ch4 + hv -> 1ch2 + h2 |
---|
340 | ! j_ch4_3ch2_h_h = 14 ! ch4 + hv -> 3ch2 + h + h |
---|
341 | ! j_ch4_ch_h2_h = 15 ! ch4 + hv -> ch + h2 + h |
---|
342 | ! j_ch3o2h = 16 ! ch3o2h + hv -> ch3o + oh |
---|
343 | ! j_ch2o_hco = 17 ! ch2o + hv -> h + hco |
---|
344 | ! j_ch2o_co = 18 ! ch2o + hv -> h2 + co |
---|
345 | ! j_ch3oh = 19 ! ch3oh + hv -> ch3o + h |
---|
346 | ! j_c2h6 = 20 ! c2h6 + hv -> products |
---|
347 | ! j_hcl = 21 ! hcl + hv -> h + cl |
---|
348 | ! j_hocl = 22 ! hocl + hv -> oh + cl |
---|
349 | ! j_clo = 23 ! clo + hv -> o + cl |
---|
350 | ! j_so2 = 24 ! so2 + hv -> so + o |
---|
351 | ! j_so = 25 ! so + hv -> s + o |
---|
352 | ! j_h2s = 26 ! h2s + hv -> hs + s |
---|
353 | ! j_so3 = 27 ! so2 + hv -> so2 + o |
---|
354 | ! j_hno3 = 28 ! hno3 + hv -> oh + no2 |
---|
355 | ! j_hno4 = 29 ! hno4 + hv -> ho2 + no2 |
---|
356 | |
---|
357 | ! fill v_phot array |
---|
358 | |
---|
359 | v_phot(:,:) = 0. |
---|
360 | |
---|
361 | do l = 1,lswitch-1 |
---|
362 | v_phot(l, 1) = j(l,j_o2_o) |
---|
363 | v_phot(l, 2) = j(l,j_o2_o1d) |
---|
364 | v_phot(l, 3) = j(l,j_co2_o) |
---|
365 | v_phot(l, 4) = j(l,j_co2_o1d) |
---|
366 | v_phot(l, 5) = j(l,j_o3_o1d) |
---|
367 | v_phot(l, 6) = j(l,j_o3_o) |
---|
368 | v_phot(l, 7) = j(l,j_h2o) |
---|
369 | v_phot(l, 8) = j(l,j_h2o2) |
---|
370 | v_phot(l, 9) = j(l,j_ho2) |
---|
371 | v_phot(l,10) = j(l,j_no) |
---|
372 | v_phot(l,11) = j(l,j_no2) |
---|
373 | end do |
---|
374 | |
---|
375 | return |
---|
376 | end |
---|
377 | |
---|
378 | !***************************************************************** |
---|