1 | !============================================================================== |
---|
2 | |
---|
3 | subroutine photolysis_online(nlayer, deutchem, nb_phot_max, |
---|
4 | $ alt, press, temp, zmmean, |
---|
5 | $ i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, |
---|
6 | $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, |
---|
7 | $ i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, |
---|
8 | $ tau, sza, dist_sol, v_phot) |
---|
9 | |
---|
10 | use photolysis_mod |
---|
11 | |
---|
12 | implicit none |
---|
13 | |
---|
14 | ! input |
---|
15 | |
---|
16 | logical, intent(in) :: deutchem |
---|
17 | integer, intent(in) :: nesp ! total number of chemical species |
---|
18 | integer, intent(in) :: nlayer |
---|
19 | integer, intent(in) :: nb_phot_max |
---|
20 | integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, |
---|
21 | $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, |
---|
22 | $ i_n, i_n2d, i_no, i_no2, i_n2 |
---|
23 | |
---|
24 | real, dimension(nlayer), intent(in) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1) |
---|
25 | real, dimension(nlayer), intent(in) :: alt ! altitude (km) |
---|
26 | real, dimension(nlayer,nesp), intent(in) :: rm ! mixing ratios |
---|
27 | real, intent(in) :: tau ! integrated aerosol optical depth at the surface |
---|
28 | real, intent(in) :: sza ! solar zenith angle (degrees) |
---|
29 | real, intent(in) :: dist_sol ! solar distance (au) |
---|
30 | |
---|
31 | ! output |
---|
32 | |
---|
33 | real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1) |
---|
34 | |
---|
35 | ! solar flux at mars |
---|
36 | |
---|
37 | real, dimension(nw) :: fmars ! solar flux (w.m-2.nm-1) |
---|
38 | real :: factor |
---|
39 | |
---|
40 | ! cross-sections |
---|
41 | |
---|
42 | real, dimension(nlayer,nw,nphot) :: sj ! general cross-section array (cm2) |
---|
43 | |
---|
44 | ! atmosphere |
---|
45 | |
---|
46 | real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) |
---|
47 | real, dimension(nlayer,nw) :: dtrl ! rayleigh optical depth |
---|
48 | real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth |
---|
49 | real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo |
---|
50 | real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter |
---|
51 | real, dimension(nlayer,nw) :: dtcld ! cloud optical depth |
---|
52 | real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo |
---|
53 | real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter |
---|
54 | real, dimension(nlayer,nw,nabs) :: dtgas ! optical depth for each gas |
---|
55 | real, dimension(nlayer,nw) :: dagas ! total gas optical depth |
---|
56 | real, dimension(nlayer) :: edir, edn, eup ! normalised irradiances |
---|
57 | real, dimension(nlayer) :: fdir, fdn, fup ! normalised actinic fluxes |
---|
58 | real, dimension(nlayer) :: saflux ! total actinic flux |
---|
59 | |
---|
60 | integer, dimension(0:nlayer) :: nid |
---|
61 | real, dimension(0:nlayer,nlayer) :: dsdh |
---|
62 | |
---|
63 | integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o, |
---|
64 | $ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2, |
---|
65 | $ j_hdo_od, j_hdo_d |
---|
66 | |
---|
67 | integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no, |
---|
68 | $ a_no2, a_n2 |
---|
69 | integer :: i, ilay, iw, ialt |
---|
70 | real :: deltaj |
---|
71 | |
---|
72 | ! absorbing gas numbering |
---|
73 | |
---|
74 | a_o2 = 1 ! o2 |
---|
75 | a_co2 = 2 ! co2 |
---|
76 | a_o3 = 3 ! o3 |
---|
77 | a_h2o = 4 ! h2o |
---|
78 | a_h2o2 = 5 ! h2o2 |
---|
79 | a_ho2 = 6 ! ho2 |
---|
80 | a_h2 = 7 ! h2 |
---|
81 | a_no = 8 ! no |
---|
82 | a_no2 = 9 ! no2 |
---|
83 | a_n2 = 10 ! n2 |
---|
84 | |
---|
85 | ! photodissociation rates numbering. |
---|
86 | ! photodissociations must be ordered the same way in subroutine "indice" |
---|
87 | |
---|
88 | j_o2_o = 1 ! o2 + hv -> o + o |
---|
89 | j_o2_o1d = 2 ! o2 + hv -> o + o(1d) |
---|
90 | j_co2_o = 3 ! co2 + hv -> co + o |
---|
91 | j_co2_o1d = 4 ! co2 + hv -> co + o(1d) |
---|
92 | j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d) |
---|
93 | j_o3_o = 6 ! o3 + hv -> o2 + o |
---|
94 | j_h2o = 7 ! h2o + hv -> h + oh |
---|
95 | j_h2o2 = 8 ! h2o2 + hv -> oh + oh |
---|
96 | j_ho2 = 9 ! ho2 + hv -> oh + o |
---|
97 | j_h2 = 10 ! h2 + hv -> h + h |
---|
98 | j_no = 11 ! no + hv -> n + o |
---|
99 | j_no2 = 12 ! no2 + hv -> no + o |
---|
100 | j_n2 = 13 ! n2 + hv -> n + n |
---|
101 | j_hdo_od = 14 ! hdo + hv -> od + h |
---|
102 | j_hdo_d = 15 ! hdo + hv -> d + oh |
---|
103 | |
---|
104 | !==== air column increments and rayleigh optical depth |
---|
105 | |
---|
106 | call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl) |
---|
107 | |
---|
108 | !==== set temperature-dependent cross-sections and optical depths |
---|
109 | |
---|
110 | ! o2: |
---|
111 | |
---|
112 | call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200, |
---|
113 | $ xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d, |
---|
114 | $ colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj) |
---|
115 | |
---|
116 | ! co2: |
---|
117 | |
---|
118 | call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295, |
---|
119 | $ xsco2_370, yieldco2, j_co2_o, j_co2_o1d, |
---|
120 | $ colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj) |
---|
121 | |
---|
122 | ! o3: |
---|
123 | |
---|
124 | call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298, |
---|
125 | $ j_o3_o, j_o3_o1d, |
---|
126 | $ colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj) |
---|
127 | |
---|
128 | ! h2o2: |
---|
129 | |
---|
130 | call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2, |
---|
131 | $ colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj) |
---|
132 | |
---|
133 | ! no2: |
---|
134 | |
---|
135 | call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220, |
---|
136 | $ xsno2_294, yldno2_248, yldno2_298, j_no2, |
---|
137 | $ colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj) |
---|
138 | |
---|
139 | !==== temperature independent optical depths and cross-sections |
---|
140 | |
---|
141 | ! optical depths |
---|
142 | |
---|
143 | do ilay = 1,nlayer |
---|
144 | do iw = 1,nw-1 |
---|
145 | dtgas(ilay,iw,a_h2o) = colinc(ilay)*rm(ilay,i_h2o)*xsh2o(iw) |
---|
146 | dtgas(ilay,iw,a_ho2) = colinc(ilay)*rm(ilay,i_ho2)*xsho2(iw) |
---|
147 | dtgas(ilay,iw,a_h2) = colinc(ilay)*rm(ilay,i_h2)*xsh2(iw) |
---|
148 | dtgas(ilay,iw,a_no) = colinc(ilay)*rm(ilay,i_no)*xsno(iw) |
---|
149 | dtgas(ilay,iw,a_n2) = colinc(ilay)*rm(ilay,i_n2)*xsn2(iw) |
---|
150 | end do |
---|
151 | end do |
---|
152 | |
---|
153 | ! total gas optical depth |
---|
154 | |
---|
155 | dagas(:,:) = 0. |
---|
156 | |
---|
157 | do ilay = 1,nlayer |
---|
158 | do iw = 1,nw-1 |
---|
159 | do i = 1,nabs |
---|
160 | dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i) |
---|
161 | end do |
---|
162 | end do |
---|
163 | end do |
---|
164 | |
---|
165 | ! cross-sections |
---|
166 | |
---|
167 | do ilay = 1,nlayer |
---|
168 | do iw = 1,nw-1 |
---|
169 | sj(ilay,iw,j_h2o) = xsh2o(iw) ! h2o |
---|
170 | sj(ilay,iw,j_ho2) = xsho2(iw) ! ho2 |
---|
171 | sj(ilay,iw,j_h2) = xsh2(iw)*yieldh2(iw) ! h2 |
---|
172 | sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no |
---|
173 | sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2 |
---|
174 | end do |
---|
175 | end do |
---|
176 | |
---|
177 | !HDO cross section |
---|
178 | if (deutchem) then |
---|
179 | do ilay=1,nlayer |
---|
180 | do iw=1,nw-1 |
---|
181 | !Two chanels for HDO, with same cross section |
---|
182 | sj(ilay,iw,j_hdo_od) = 0.5*xshdo(iw) ! hdo -> od + h |
---|
183 | sj(ilay,iw,j_hdo_d) = 0.5*xshdo(iw) ! hdo -> d + oh |
---|
184 | enddo |
---|
185 | enddo |
---|
186 | endif |
---|
187 | |
---|
188 | !==== set aerosol properties and optical depth |
---|
189 | |
---|
190 | call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer) |
---|
191 | |
---|
192 | !==== set cloud properties and optical depth |
---|
193 | |
---|
194 | call setcld(nlayer,nw,dtcld,omcld,gcld) |
---|
195 | |
---|
196 | !==== slant path lengths in spherical geometry |
---|
197 | |
---|
198 | call sphers(nlayer,alt,sza,dsdh,nid) |
---|
199 | |
---|
200 | !==== solar flux at mars |
---|
201 | |
---|
202 | factor = (1./dist_sol)**2. |
---|
203 | do iw = 1,nw-1 |
---|
204 | fmars(iw) = f(iw)*factor |
---|
205 | end do |
---|
206 | |
---|
207 | !==== initialise photolysis rates |
---|
208 | |
---|
209 | v_phot(:,1:nphot) = 0. |
---|
210 | |
---|
211 | !==== start of wavelength lopp |
---|
212 | |
---|
213 | do iw = 1,nw-1 |
---|
214 | |
---|
215 | ! monochromatic radiative transfer. outputs are: |
---|
216 | ! normalized irradiances edir(nlayer), edn(nlayer), eup(nlayer) |
---|
217 | ! normalized actinic fluxes fdir(nlayer), fdn(nlayer), fup(nlayer) |
---|
218 | ! where |
---|
219 | ! dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse |
---|
220 | |
---|
221 | call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl, |
---|
222 | $ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer, |
---|
223 | $ edir, edn, eup, fdir, fdn, fup) |
---|
224 | |
---|
225 | |
---|
226 | ! spherical actinic flux |
---|
227 | |
---|
228 | do ilay = 1,nlayer |
---|
229 | saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay)) |
---|
230 | end do |
---|
231 | |
---|
232 | ! photolysis rate integration |
---|
233 | |
---|
234 | do i = 1,nphot |
---|
235 | do ilay = 1,nlayer |
---|
236 | deltaj = saflux(ilay)*sj(ilay,iw,i) |
---|
237 | v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw)) |
---|
238 | end do |
---|
239 | end do |
---|
240 | |
---|
241 | ! eliminate small values |
---|
242 | |
---|
243 | where (v_phot(:,1:nphot) < 1.e-30) |
---|
244 | v_phot(:,1:nphot) = 0. |
---|
245 | end where |
---|
246 | |
---|
247 | end do ! iw |
---|
248 | |
---|
249 | contains |
---|
250 | |
---|
251 | !============================================================================== |
---|
252 | |
---|
253 | subroutine setair(nlev, nw, wl, wc, press, temp, zmmean, |
---|
254 | $ colinc, dtrl) |
---|
255 | |
---|
256 | *-----------------------------------------------------------------------------* |
---|
257 | *= PURPOSE: =* |
---|
258 | *= computes air column increments and rayleigh optical depth =* |
---|
259 | *-----------------------------------------------------------------------------* |
---|
260 | |
---|
261 | implicit none |
---|
262 | |
---|
263 | ! input: |
---|
264 | |
---|
265 | integer :: nlev, nw |
---|
266 | |
---|
267 | real, dimension(nw) :: wl, wc ! lower and central wavelength grid (nm) |
---|
268 | real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1) |
---|
269 | |
---|
270 | ! output: |
---|
271 | |
---|
272 | real, dimension(nlev) :: colinc ! air column increments (molecule.cm-2) |
---|
273 | real, dimension(nlev,nw) :: dtrl ! rayleigh optical depth |
---|
274 | |
---|
275 | ! local: |
---|
276 | |
---|
277 | real, parameter :: avo = 6.022e23 |
---|
278 | real, parameter :: g = 3.72 |
---|
279 | real :: dp, nu |
---|
280 | real, dimension(nw) :: srayl |
---|
281 | integer :: ilev, iw |
---|
282 | |
---|
283 | ! compute column increments |
---|
284 | |
---|
285 | do ilev = 1, nlev-1 |
---|
286 | dp = (press(ilev) - press(ilev+1))*100. |
---|
287 | colinc(ilev) = avo*0.1*dp/(zmmean(ilev)*g) |
---|
288 | end do |
---|
289 | colinc(nlev) = avo*0.1*press(nlev)*100./(zmmean(nlev)*g) |
---|
290 | |
---|
291 | do iw = 1, nw - 1 |
---|
292 | |
---|
293 | ! co2 rayleigh cross-section |
---|
294 | ! ityaksov et al., chem. phys. lett., 462, 31-34, 2008 |
---|
295 | |
---|
296 | nu = 1./(wc(iw)*1.e-7) |
---|
297 | srayl(iw) = 1.78e-26*nu**(4. + 0.625) |
---|
298 | srayl(iw) = srayl(iw)*1.e-20 ! cm2 |
---|
299 | |
---|
300 | do ilev = 1, nlev |
---|
301 | dtrl(ilev,iw) = colinc(ilev)*srayl(iw) ! cm2 |
---|
302 | end do |
---|
303 | end do |
---|
304 | |
---|
305 | end subroutine setair |
---|
306 | |
---|
307 | !============================================================================== |
---|
308 | |
---|
309 | subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295, |
---|
310 | $ xsco2_370, yieldco2, j_co2_o, j_co2_o1d, |
---|
311 | $ colinc, rm, dt, sj) |
---|
312 | |
---|
313 | !-----------------------------------------------------------------------------* |
---|
314 | != PURPOSE: =* |
---|
315 | != Set up the CO2 temperature-dependent cross-sections and optical depth =* |
---|
316 | !-----------------------------------------------------------------------------* |
---|
317 | |
---|
318 | implicit none |
---|
319 | |
---|
320 | ! input: |
---|
321 | |
---|
322 | integer :: nd ! number of photolysis rates |
---|
323 | integer :: nlayer ! number of vertical layers |
---|
324 | integer :: nw ! number of wavelength grid points |
---|
325 | integer :: j_co2_o, j_co2_o1d ! photolysis indexes |
---|
326 | real, dimension(nw) :: wc ! central wavelength for each interval |
---|
327 | real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2) |
---|
328 | real, dimension(nw) :: yieldco2 ! co2 photodissociation yield |
---|
329 | real, dimension(nlayer) :: tlay ! temperature (k) |
---|
330 | real, dimension(nlayer) :: rm ! co2 mixing ratio |
---|
331 | real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) |
---|
332 | |
---|
333 | ! output: |
---|
334 | |
---|
335 | real, dimension(nlayer,nw) :: dt ! optical depth |
---|
336 | real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) |
---|
337 | |
---|
338 | ! local: |
---|
339 | |
---|
340 | integer :: extrapol |
---|
341 | integer :: i, l |
---|
342 | real :: temp, sco2 |
---|
343 | |
---|
344 | ! extrapol = 0 no extrapolation below 195 k |
---|
345 | ! extrapol = 1 extrapolation below 195 k |
---|
346 | |
---|
347 | extrapol = 0 |
---|
348 | |
---|
349 | do i = 1, nlayer |
---|
350 | if (extrapol == 1) then |
---|
351 | temp = tlay(i) |
---|
352 | else |
---|
353 | temp = max(tlay(i), 195.) |
---|
354 | end if |
---|
355 | temp = min(temp, 370.) |
---|
356 | do l = 1, nw-1 |
---|
357 | if (temp <= 295.) then |
---|
358 | if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then |
---|
359 | sco2 = alog(xsco2_195(l)) |
---|
360 | $ + (alog(xsco2_295(l)) - alog(xsco2_195(l))) |
---|
361 | $ /(295. - 195.)*(temp - 195.) |
---|
362 | sco2 = exp(sco2) |
---|
363 | else |
---|
364 | sco2 = 0. |
---|
365 | end if |
---|
366 | else |
---|
367 | if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then |
---|
368 | sco2 = alog(xsco2_295(l)) |
---|
369 | $ + (alog(xsco2_370(l)) - alog(xsco2_295(l))) |
---|
370 | $ /(370. - 295.)*(temp - 295.) |
---|
371 | sco2 = exp(sco2) |
---|
372 | else |
---|
373 | sco2 = 0. |
---|
374 | end if |
---|
375 | end if |
---|
376 | |
---|
377 | ! optical depth |
---|
378 | |
---|
379 | dt(i,l) = colinc(i)*rm(i)*sco2 |
---|
380 | |
---|
381 | ! production of o(1d) for wavelengths shorter than 167 nm |
---|
382 | |
---|
383 | if (wc(l) >= 167.) then |
---|
384 | sj(i,l,j_co2_o) = sco2*yieldco2(l) |
---|
385 | sj(i,l,j_co2_o1d) = 0. |
---|
386 | else |
---|
387 | sj(i,l,j_co2_o) = 0. |
---|
388 | sj(i,l,j_co2_o1d) = sco2*yieldco2(l) |
---|
389 | end if |
---|
390 | end do |
---|
391 | end do |
---|
392 | |
---|
393 | end subroutine setco2 |
---|
394 | |
---|
395 | !============================================================================== |
---|
396 | |
---|
397 | subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150, |
---|
398 | $ xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o, |
---|
399 | $ j_o2_o1d, colinc, rm, dt, sj) |
---|
400 | |
---|
401 | !-----------------------------------------------------------------------------* |
---|
402 | != PURPOSE: =* |
---|
403 | != Set up the O2 temperature-dependent cross-sections and optical depth =* |
---|
404 | !-----------------------------------------------------------------------------* |
---|
405 | |
---|
406 | implicit none |
---|
407 | |
---|
408 | ! input: |
---|
409 | |
---|
410 | integer :: nd ! number of photolysis rates |
---|
411 | integer :: nlayer ! number of vertical layers |
---|
412 | integer :: nw ! number of wavelength grid points |
---|
413 | integer :: mopt ! high-res/low-res switch |
---|
414 | integer :: j_o2_o, j_o2_o1d ! photolysis indexes |
---|
415 | real, dimension(nw) :: wc ! central wavelength for each interval |
---|
416 | real, dimension(nw) :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2) |
---|
417 | $ xso2_300 |
---|
418 | real, dimension(nw) :: yieldo2 ! o2 photodissociation yield |
---|
419 | real, dimension(nlayer) :: tlay ! temperature (k) |
---|
420 | real, dimension(nlayer) :: rm ! o2 mixing ratio |
---|
421 | real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) |
---|
422 | |
---|
423 | ! output: |
---|
424 | |
---|
425 | real, dimension(nlayer,nw) :: dt ! optical depth |
---|
426 | real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) |
---|
427 | |
---|
428 | ! local: |
---|
429 | |
---|
430 | integer :: ilev, iw |
---|
431 | real :: temp |
---|
432 | real :: xso2, factor |
---|
433 | |
---|
434 | ! correction by factor if low-resolution in schumann-runge bands |
---|
435 | |
---|
436 | if (mopt == 1) then |
---|
437 | factor = 1. |
---|
438 | else if (mopt == 2) then |
---|
439 | factor = 0.8 |
---|
440 | end if |
---|
441 | |
---|
442 | ! calculate temperature dependance |
---|
443 | |
---|
444 | do ilev = 1,nlayer |
---|
445 | temp = max(tlay(ilev),150.) |
---|
446 | temp = min(temp, 300.) |
---|
447 | do iw = 1, nw-1 |
---|
448 | if (tlay(ilev) > 250.) then |
---|
449 | xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw)) |
---|
450 | $ /(300. - 250.)*(temp - 250.) |
---|
451 | else if (tlay(ilev) > 200.) then |
---|
452 | xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw)) |
---|
453 | $ /(250. - 200.)*(temp - 200.) |
---|
454 | else |
---|
455 | xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw)) |
---|
456 | $ /(200. - 150.)*(temp - 150.) |
---|
457 | end if |
---|
458 | |
---|
459 | if (wc(iw) > 180. .and. wc(iw) < 200.) then |
---|
460 | xso2 = xso2*factor |
---|
461 | end if |
---|
462 | |
---|
463 | ! optical depth |
---|
464 | |
---|
465 | dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso2 |
---|
466 | |
---|
467 | ! production of o(1d) for wavelengths shorter than 175 nm |
---|
468 | |
---|
469 | if (wc(iw) >= 175.) then |
---|
470 | sj(ilev,iw,j_o2_o) = xso2*yieldo2(iw) |
---|
471 | sj(ilev,iw,j_o2_o1d) = 0. |
---|
472 | else |
---|
473 | sj(ilev,iw,j_o2_o) = 0. |
---|
474 | sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw) |
---|
475 | end if |
---|
476 | |
---|
477 | end do |
---|
478 | end do |
---|
479 | |
---|
480 | end subroutine seto2 |
---|
481 | |
---|
482 | !============================================================================== |
---|
483 | |
---|
484 | subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298, |
---|
485 | $ j_o3_o, j_o3_o1d, |
---|
486 | $ colinc, rm, dt, sj) |
---|
487 | |
---|
488 | !-----------------------------------------------------------------------------* |
---|
489 | != PURPOSE: =* |
---|
490 | != Set up the O3 temperature dependent cross-sections and optical depth =* |
---|
491 | !-----------------------------------------------------------------------------* |
---|
492 | |
---|
493 | implicit none |
---|
494 | |
---|
495 | ! input: |
---|
496 | |
---|
497 | integer :: nd ! number of photolysis rates |
---|
498 | integer :: nlayer ! number of vertical layers |
---|
499 | integer :: nw ! number of wavelength grid points |
---|
500 | integer :: j_o3_o, j_o3_o1d ! photolysis indexes |
---|
501 | real, dimension(nw) :: wc ! central wavelength for each interval |
---|
502 | real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2) |
---|
503 | real, dimension(nlayer) :: tlay ! temperature (k) |
---|
504 | real, dimension(nlayer) :: rm ! o3 mixing ratio |
---|
505 | real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) |
---|
506 | |
---|
507 | ! output: |
---|
508 | |
---|
509 | real, dimension(nlayer,nw) :: dt ! optical depth |
---|
510 | real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) |
---|
511 | |
---|
512 | ! local: |
---|
513 | ! |
---|
514 | integer :: ilev, iw |
---|
515 | real :: temp |
---|
516 | real, dimension(nw) :: xso3(nw) |
---|
517 | real, dimension(nw) :: qy1d ! quantum yield for o(1d) production |
---|
518 | real :: q1, q2, a1, a2, a3 |
---|
519 | |
---|
520 | do ilev = 1, nlayer |
---|
521 | temp = max(tlay(ilev), 218.) |
---|
522 | temp = min(temp,298.) |
---|
523 | do iw = 1, nw-1 |
---|
524 | xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw)) |
---|
525 | $ /(298. - 218.) *(temp - 218.) |
---|
526 | |
---|
527 | ! optical depth |
---|
528 | |
---|
529 | dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso3(iw) |
---|
530 | |
---|
531 | end do |
---|
532 | |
---|
533 | ! calculate quantum yield for o(1d) production (jpl 2006) |
---|
534 | |
---|
535 | temp = max(tlay(ilev),200.) |
---|
536 | temp = min(temp,320.) |
---|
537 | do iw = 1, nw-1 |
---|
538 | if (wc(iw) <= 306.) then |
---|
539 | qy1d(iw) = 0.90 |
---|
540 | else if (wc(iw) > 306. .and. wc(iw) < 328.) then |
---|
541 | q1 = 1. |
---|
542 | q2 = exp(-825.518/(0.695*temp)) |
---|
543 | a1 = (304.225 - wc(iw))/5.576 |
---|
544 | a2 = (314.957 - wc(iw))/6.601 |
---|
545 | a3 = (310.737 - wc(iw))/2.187 |
---|
546 | qy1d(iw) = (q1/(q1 + q2))*0.8036*exp(-(a1*a1*a1*a1)) |
---|
547 | $ + (q2/(q1 + q2))*8.9061*(temp/300.)**2. |
---|
548 | $ *exp(-(a2*a2)) |
---|
549 | $ + 0.1192*(temp/300.)**1.5*exp(-(a3*a3)) |
---|
550 | $ + 0.0765 |
---|
551 | else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then |
---|
552 | qy1d(iw) = 0.08 |
---|
553 | else |
---|
554 | qy1d(iw) = 0. |
---|
555 | endif |
---|
556 | end do |
---|
557 | do iw = 1, nw-1 |
---|
558 | sj(ilev,iw,j_o3_o) = xso3(iw)*(1. - qy1d(iw)) |
---|
559 | sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw) |
---|
560 | end do |
---|
561 | end do |
---|
562 | |
---|
563 | end subroutine seto3 |
---|
564 | |
---|
565 | !============================================================================== |
---|
566 | |
---|
567 | subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2, |
---|
568 | $ colinc, rm, dt, sj) |
---|
569 | |
---|
570 | !-----------------------------------------------------------------------------* |
---|
571 | != PURPOSE: =* |
---|
572 | != Set up the h2o2 temperature dependent cross-sections and optical depth =* |
---|
573 | !-----------------------------------------------------------------------------* |
---|
574 | |
---|
575 | implicit none |
---|
576 | |
---|
577 | ! input: |
---|
578 | |
---|
579 | integer :: nd ! number of photolysis rates |
---|
580 | integer :: nlayer ! number of vertical layers |
---|
581 | integer :: nw ! number of wavelength grid points |
---|
582 | integer :: j_h2o2 ! photolysis index |
---|
583 | real, dimension(nw) :: wc ! central wavelength for each interval |
---|
584 | real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2) |
---|
585 | real, dimension(nlayer) :: tlay ! temperature (k) |
---|
586 | real, dimension(nlayer) :: rm ! h2o2 mixing ratio |
---|
587 | real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) |
---|
588 | |
---|
589 | ! output: |
---|
590 | |
---|
591 | real, dimension(nlayer,nw) :: dt ! optical depth |
---|
592 | real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) |
---|
593 | |
---|
594 | ! local: |
---|
595 | |
---|
596 | integer :: ilev, iw |
---|
597 | real :: a0, a1, a2, a3, a4, a5, a6, a7 |
---|
598 | real :: b0, b1, b2, b3, b4 |
---|
599 | real :: lambda, suma, sumb, chi, temp, xs |
---|
600 | |
---|
601 | A0 = 6.4761E+04 |
---|
602 | A1 = -9.2170972E+02 |
---|
603 | A2 = 4.535649 |
---|
604 | A3 = -4.4589016E-03 |
---|
605 | A4 = -4.035101E-05 |
---|
606 | A5 = 1.6878206E-07 |
---|
607 | A6 = -2.652014E-10 |
---|
608 | A7 = 1.5534675E-13 |
---|
609 | |
---|
610 | B0 = 6.8123E+03 |
---|
611 | B1 = -5.1351E+01 |
---|
612 | B2 = 1.1522E-01 |
---|
613 | B3 = -3.0493E-05 |
---|
614 | B4 = -1.0924E-07 |
---|
615 | |
---|
616 | ! temperature dependance: jpl 2006 |
---|
617 | |
---|
618 | do ilev = 1,nlayer |
---|
619 | temp = min(max(tlay(ilev),200.),400.) |
---|
620 | chi = 1./(1. + exp(-1265./temp)) |
---|
621 | do iw = 1, nw-1 |
---|
622 | if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then |
---|
623 | lambda = wc(iw) |
---|
624 | sumA = ((((((A7*lambda + A6)*lambda + A5)*lambda + |
---|
625 | $ A4)*lambda +A3)*lambda + A2)*lambda + |
---|
626 | $ A1)*lambda + A0 |
---|
627 | sumB = (((B4*lambda + B3)*lambda + B2)*lambda + |
---|
628 | $ B1)*lambda + B0 |
---|
629 | xs = (chi*sumA + (1. - chi)*sumB)*1.e-21 |
---|
630 | sj(ilev,iw,j_h2o2) = xs |
---|
631 | else |
---|
632 | sj(ilev,iw,j_h2o2) = xsh2o2(iw) |
---|
633 | end if |
---|
634 | |
---|
635 | ! optical depth |
---|
636 | |
---|
637 | dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_h2o2) |
---|
638 | end do |
---|
639 | end do |
---|
640 | |
---|
641 | end subroutine seth2o2 |
---|
642 | |
---|
643 | !============================================================================== |
---|
644 | |
---|
645 | subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220, |
---|
646 | $ xsno2_294, yldno2_248, yldno2_298, j_no2, |
---|
647 | $ colinc, rm, dt, sj) |
---|
648 | |
---|
649 | !-----------------------------------------------------------------------------* |
---|
650 | != PURPOSE: =* |
---|
651 | != Set up the no2 temperature-dependent cross-sections and optical depth =* |
---|
652 | !-----------------------------------------------------------------------------* |
---|
653 | |
---|
654 | implicit none |
---|
655 | |
---|
656 | ! input: |
---|
657 | |
---|
658 | integer :: nd ! number of photolysis rates |
---|
659 | integer :: nlayer ! number of vertical layers |
---|
660 | integer :: nw ! number of wavelength grid points |
---|
661 | integer :: j_no2 ! photolysis index |
---|
662 | real, dimension(nw) :: wc ! central wavelength for each interval |
---|
663 | real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2) |
---|
664 | real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k |
---|
665 | real, dimension(nlayer) :: tlay ! temperature (k) |
---|
666 | real, dimension(nlayer) :: rm ! no2 mixing ratio |
---|
667 | real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2) |
---|
668 | |
---|
669 | ! output: |
---|
670 | |
---|
671 | real, dimension(nlayer,nw) :: dt ! optical depth |
---|
672 | real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2) |
---|
673 | |
---|
674 | ! local: |
---|
675 | |
---|
676 | integer :: ilev, iw |
---|
677 | real :: temp, qy |
---|
678 | |
---|
679 | ! temperature dependance: jpl 2006 |
---|
680 | |
---|
681 | do ilev = 1,nlayer |
---|
682 | temp = max(220.,min(tlay(ilev),294.)) |
---|
683 | do iw = 1, nw - 1 |
---|
684 | if (wc(iw) < 238.) then |
---|
685 | sj(ilev,iw,j_no2) = xsno2(iw) |
---|
686 | else |
---|
687 | sj(ilev,iw,j_no2) = xsno2_220(iw) |
---|
688 | $ + (xsno2_294(iw) - xsno2_220(iw)) |
---|
689 | $ /(294. - 220.)*(temp - 220.) |
---|
690 | end if |
---|
691 | |
---|
692 | ! optical depth |
---|
693 | |
---|
694 | dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_no2) |
---|
695 | end do |
---|
696 | end do |
---|
697 | |
---|
698 | ! quantum yield: jpl 2006 |
---|
699 | |
---|
700 | do ilev = 1,nlayer |
---|
701 | temp = max(248.,min(tlay(ilev),298.)) |
---|
702 | do iw = 1, nw - 1 |
---|
703 | qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw)) |
---|
704 | $ /(298. - 248.)*(temp - 248.) |
---|
705 | sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy |
---|
706 | end do |
---|
707 | end do |
---|
708 | |
---|
709 | end subroutine setno2 |
---|
710 | |
---|
711 | !============================================================================== |
---|
712 | |
---|
713 | subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer) |
---|
714 | |
---|
715 | !-----------------------------------------------------------------------------* |
---|
716 | != PURPOSE: =* |
---|
717 | != Set aerosol properties for each specified altitude layer. Properties =* |
---|
718 | != may be wavelength dependent. =* |
---|
719 | !-----------------------------------------------------------------------------* |
---|
720 | |
---|
721 | implicit none |
---|
722 | |
---|
723 | ! input |
---|
724 | |
---|
725 | integer :: nlayer ! number of vertical layers |
---|
726 | integer :: nw ! number of wavelength grid points |
---|
727 | real, dimension(nlayer) :: alt ! altitude (km) |
---|
728 | real :: tau ! integrated aerosol optical depth at the surface |
---|
729 | |
---|
730 | ! output |
---|
731 | |
---|
732 | real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth |
---|
733 | real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo |
---|
734 | real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter |
---|
735 | |
---|
736 | ! local |
---|
737 | |
---|
738 | integer :: ilay, iw |
---|
739 | real, dimension(nlayer) :: aer ! dust extinction |
---|
740 | real :: omega, g, scaleh, gamma |
---|
741 | real :: dz, tautot, q0 |
---|
742 | |
---|
743 | omega = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm |
---|
744 | g = 0.88 ! asymmetry factor : mateshvili et al. (2007) at 210 nm |
---|
745 | scaleh = 10. ! scale height (km) |
---|
746 | gamma = 0.03 ! conrath parameter |
---|
747 | |
---|
748 | dtaer(:,:) = 0. |
---|
749 | omaer(:,:) = 0. |
---|
750 | gaer(:,:) = 0. |
---|
751 | |
---|
752 | ! optical depth profile: |
---|
753 | |
---|
754 | tautot = 0. |
---|
755 | do ilay = 1, nlayer-1 |
---|
756 | dz = alt(ilay+1) - alt(ilay) |
---|
757 | tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz |
---|
758 | end do |
---|
759 | |
---|
760 | q0 = tau/tautot |
---|
761 | do ilay = 1, nlayer-1 |
---|
762 | dz = alt(ilay+1) - alt(ilay) |
---|
763 | dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz |
---|
764 | omaer(ilay,:) = omega |
---|
765 | gaer(ilay,:) = g |
---|
766 | end do |
---|
767 | |
---|
768 | end subroutine setaer |
---|
769 | |
---|
770 | !============================================================================== |
---|
771 | |
---|
772 | subroutine setcld(nlayer,nw,dtcld,omcld,gcld) |
---|
773 | |
---|
774 | !-----------------------------------------------------------------------------* |
---|
775 | != PURPOSE: =* |
---|
776 | != Set cloud properties for each specified altitude layer. Properties =* |
---|
777 | != may be wavelength dependent. =* |
---|
778 | !-----------------------------------------------------------------------------* |
---|
779 | |
---|
780 | implicit none |
---|
781 | |
---|
782 | ! input |
---|
783 | |
---|
784 | integer :: nlayer ! number of vertical layers |
---|
785 | integer :: nw ! number of wavelength grid points |
---|
786 | |
---|
787 | ! output |
---|
788 | |
---|
789 | real, dimension(nlayer,nw) :: dtcld ! cloud optical depth |
---|
790 | real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo |
---|
791 | real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter |
---|
792 | |
---|
793 | ! local |
---|
794 | |
---|
795 | integer :: ilay, iw |
---|
796 | |
---|
797 | ! dtcld : optical depth |
---|
798 | ! omcld : single scattering albedo |
---|
799 | ! gcld : asymmetry factor |
---|
800 | |
---|
801 | do ilay = 1, nlayer - 1 |
---|
802 | do iw = 1, nw - 1 |
---|
803 | dtcld(ilay,iw) = 0. ! no clouds for the moment |
---|
804 | omcld(ilay,iw) = 0.99 |
---|
805 | gcld(ilay,iw) = 0.85 |
---|
806 | end do |
---|
807 | end do |
---|
808 | |
---|
809 | end subroutine setcld |
---|
810 | |
---|
811 | !============================================================================== |
---|
812 | |
---|
813 | subroutine sphers(nlev, z, zen, dsdh, nid) |
---|
814 | |
---|
815 | !-----------------------------------------------------------------------------* |
---|
816 | != PURPOSE: =* |
---|
817 | != Calculate slant path over vertical depth ds/dh in spherical geometry. =* |
---|
818 | != Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =* |
---|
819 | != for computing the radiation field available for photolysis and heating =* |
---|
820 | != at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =* |
---|
821 | !-----------------------------------------------------------------------------* |
---|
822 | != PARAMETERS: =* |
---|
823 | != NZ - INTEGER, number of specified altitude levels in the working (I)=* |
---|
824 | != grid =* |
---|
825 | != Z - REAL, specified altitude working grid (km) (I)=* |
---|
826 | != ZEN - REAL, solar zenith angle (degrees) (I)=* |
---|
827 | != DSDH - REAL, slant path of direct beam through each layer crossed (O)=* |
---|
828 | != when travelling from the top of the atmosphere to layer i; =* |
---|
829 | != DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* |
---|
830 | != NID - INTEGER, number of layers crossed by the direct beam when (O)=* |
---|
831 | != travelling from the top of the atmosphere to layer i; =* |
---|
832 | != NID(i), i = 0..NZ-1 =* |
---|
833 | !-----------------------------------------------------------------------------* |
---|
834 | |
---|
835 | implicit none |
---|
836 | |
---|
837 | ! input |
---|
838 | |
---|
839 | integer, intent(in) :: nlev |
---|
840 | real, dimension(nlev), intent(in) :: z |
---|
841 | real, intent(in) :: zen |
---|
842 | |
---|
843 | ! output |
---|
844 | |
---|
845 | INTEGER nid(0:nlev) |
---|
846 | REAL dsdh(0:nlev,nlev) |
---|
847 | |
---|
848 | ! more program constants |
---|
849 | |
---|
850 | REAL re, ze(nlev) |
---|
851 | REAL dr |
---|
852 | real radius |
---|
853 | parameter (radius = 3393.) |
---|
854 | |
---|
855 | ! local |
---|
856 | |
---|
857 | real :: pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm |
---|
858 | integer :: i, j, k, id, nlay |
---|
859 | |
---|
860 | REAL zd(0:nlev-1) |
---|
861 | |
---|
862 | !----------------------------------------------------------------------------- |
---|
863 | |
---|
864 | pi = acos(-1.0) |
---|
865 | dr = pi/180. |
---|
866 | zenrad = zen*dr |
---|
867 | |
---|
868 | ! number of layers: |
---|
869 | |
---|
870 | nlay = nlev - 1 |
---|
871 | |
---|
872 | ! include the elevation above sea level to the radius of Mars: |
---|
873 | |
---|
874 | re = radius + z(1) |
---|
875 | |
---|
876 | ! correspondingly z changed to the elevation above Mars surface: |
---|
877 | |
---|
878 | DO k = 1, nlev |
---|
879 | ze(k) = z(k) - z(1) |
---|
880 | END DO |
---|
881 | |
---|
882 | ! inverse coordinate of z |
---|
883 | |
---|
884 | zd(0) = ze(nlev) |
---|
885 | DO k = 1, nlay |
---|
886 | zd(k) = ze(nlev - k) |
---|
887 | END DO |
---|
888 | |
---|
889 | ! initialise dsdh(i,j), nid(i) |
---|
890 | |
---|
891 | nid(:) = 0. |
---|
892 | dsdh(:,:) = 0. |
---|
893 | |
---|
894 | ! calculate ds/dh of every layer |
---|
895 | |
---|
896 | do i = 0,nlay |
---|
897 | rpsinz = (re + zd(i))*sin(zenrad) |
---|
898 | |
---|
899 | IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN |
---|
900 | nid(i) = -1 |
---|
901 | ELSE |
---|
902 | |
---|
903 | ! Find index of layer in which the screening height lies |
---|
904 | |
---|
905 | id = i |
---|
906 | if (zen > 90.) then |
---|
907 | do j = 1,nlay |
---|
908 | IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND. |
---|
909 | $ (rpsinz .GE. ( zd(j) + re )) ) id = j |
---|
910 | end do |
---|
911 | end if |
---|
912 | |
---|
913 | do j = 1,id |
---|
914 | sm = 1.0 |
---|
915 | IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0) |
---|
916 | $ sm = -1.0 |
---|
917 | |
---|
918 | rj = re + zd(j-1) |
---|
919 | rjp1 = re + zd(j) |
---|
920 | |
---|
921 | dhj = zd(j-1) - zd(j) |
---|
922 | |
---|
923 | ga = rj*rj - rpsinz*rpsinz |
---|
924 | gb = rjp1*rjp1 - rpsinz*rpsinz |
---|
925 | |
---|
926 | ga = max(ga, 0.) |
---|
927 | gb = max(gb, 0.) |
---|
928 | |
---|
929 | IF (id.GT.i .AND. j.EQ.id) THEN |
---|
930 | dsj = sqrt(ga) |
---|
931 | ELSE |
---|
932 | dsj = sqrt(ga) - sm*sqrt(gb) |
---|
933 | END IF |
---|
934 | dsdh(i,j) = dsj/dhj |
---|
935 | end do |
---|
936 | nid(i) = id |
---|
937 | end if |
---|
938 | end do ! i = 0,nlay |
---|
939 | |
---|
940 | end subroutine sphers |
---|
941 | |
---|
942 | !============================================================================== |
---|
943 | |
---|
944 | SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl, |
---|
945 | $ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer, |
---|
946 | $ edir, edn, eup, fdir, fdn, fup) |
---|
947 | |
---|
948 | implicit none |
---|
949 | |
---|
950 | ! input |
---|
951 | |
---|
952 | integer, intent(in) :: nlev, nw, iw ! number of wavelength grid points |
---|
953 | REAL ag |
---|
954 | REAL zen |
---|
955 | REAL dsdh(0:nlev,nlev) |
---|
956 | INTEGER nid(0:nlev) |
---|
957 | |
---|
958 | REAL dtrl(nlev,nw) |
---|
959 | REAL dagas(nlev,nw) |
---|
960 | REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw) |
---|
961 | REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw) |
---|
962 | |
---|
963 | ! output |
---|
964 | |
---|
965 | REAL edir(nlev), edn(nlev), eup(nlev) |
---|
966 | REAL fdir(nlev), fdn(nlev), fup(nlev) |
---|
967 | |
---|
968 | ! local: |
---|
969 | |
---|
970 | REAL dt(nlev), om(nlev), g(nlev) |
---|
971 | REAL dtabs,dtsct,dscld,dsaer,dacld,daaer |
---|
972 | INTEGER i, ii |
---|
973 | real, parameter :: largest = 1.e+36 |
---|
974 | |
---|
975 | ! specific two ps2str |
---|
976 | |
---|
977 | REAL ediri(nlev), edni(nlev), eupi(nlev) |
---|
978 | REAL fdiri(nlev), fdni(nlev), fupi(nlev) |
---|
979 | |
---|
980 | logical, save :: delta = .true. |
---|
981 | |
---|
982 | !_______________________________________________________________________ |
---|
983 | |
---|
984 | ! initialize: |
---|
985 | |
---|
986 | do i = 1, nlev |
---|
987 | fdir(i) = 0. |
---|
988 | fup(i) = 0. |
---|
989 | fdn(i) = 0. |
---|
990 | edir(i) = 0. |
---|
991 | eup(i) = 0. |
---|
992 | edn(i) = 0. |
---|
993 | end do |
---|
994 | |
---|
995 | do i = 1, nlev - 1 |
---|
996 | dscld = dtcld(i,iw)*omcld(i,iw) |
---|
997 | dacld = dtcld(i,iw)*(1.-omcld(i,iw)) |
---|
998 | |
---|
999 | dsaer = dtaer(i,iw)*omaer(i,iw) |
---|
1000 | daaer = dtaer(i,iw)*(1.-omaer(i,iw)) |
---|
1001 | |
---|
1002 | dtsct = dtrl(i,iw) + dscld + dsaer |
---|
1003 | dtabs = dagas(i,iw) + dacld + daaer |
---|
1004 | |
---|
1005 | dtabs = amax1(dtabs,1./largest) |
---|
1006 | dtsct = amax1(dtsct,1./largest) |
---|
1007 | |
---|
1008 | ! invert z-coordinate: |
---|
1009 | |
---|
1010 | ii = nlev - i |
---|
1011 | dt(ii) = dtsct + dtabs |
---|
1012 | om(ii) = dtsct/(dtsct + dtabs) |
---|
1013 | IF(dtsct .EQ. 1./largest) om(ii) = 1./largest |
---|
1014 | g(ii) = (gcld(i,iw)*dscld + |
---|
1015 | $ gaer(i,iw)*dsaer)/dtsct |
---|
1016 | end do |
---|
1017 | |
---|
1018 | ! call rt routine: |
---|
1019 | |
---|
1020 | call ps2str(nlev, zen, ag, dt, om, g, |
---|
1021 | $ dsdh, nid, delta, |
---|
1022 | $ fdiri, fupi, fdni, ediri, eupi, edni) |
---|
1023 | |
---|
1024 | ! output (invert z-coordinate) |
---|
1025 | |
---|
1026 | do i = 1, nlev |
---|
1027 | ii = nlev - i + 1 |
---|
1028 | fdir(i) = fdiri(ii) |
---|
1029 | fup(i) = fupi(ii) |
---|
1030 | fdn(i) = fdni(ii) |
---|
1031 | edir(i) = ediri(ii) |
---|
1032 | eup(i) = eupi(ii) |
---|
1033 | edn(i) = edni(ii) |
---|
1034 | end do |
---|
1035 | |
---|
1036 | end subroutine rtlink |
---|
1037 | |
---|
1038 | *=============================================================================* |
---|
1039 | |
---|
1040 | subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu, |
---|
1041 | $ dsdh, nid, delta, |
---|
1042 | $ fdr, fup, fdn, edr, eup, edn) |
---|
1043 | |
---|
1044 | !-----------------------------------------------------------------------------* |
---|
1045 | != PURPOSE: =* |
---|
1046 | != Solve two-stream equations for multiple layers. The subroutine is based =* |
---|
1047 | != on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=* |
---|
1048 | != It contains 9 two-stream methods to choose from. A pseudo-spherical =* |
---|
1049 | != correction has also been added. =* |
---|
1050 | !-----------------------------------------------------------------------------* |
---|
1051 | != PARAMETERS: =* |
---|
1052 | != NLEVEL - INTEGER, number of specified altitude levels in the working (I)=* |
---|
1053 | != grid =* |
---|
1054 | != ZEN - REAL, solar zenith angle (degrees) (I)=* |
---|
1055 | != RSFC - REAL, surface albedo at current wavelength (I)=* |
---|
1056 | != TAUU - REAL, unscaled optical depth of each layer (I)=* |
---|
1057 | != OMU - REAL, unscaled single scattering albedo of each layer (I)=* |
---|
1058 | != GU - REAL, unscaled asymmetry parameter of each layer (I)=* |
---|
1059 | != DSDH - REAL, slant path of direct beam through each layer crossed (I)=* |
---|
1060 | != when travelling from the top of the atmosphere to layer i; =* |
---|
1061 | != DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =* |
---|
1062 | != NID - INTEGER, number of layers crossed by the direct beam when (I)=* |
---|
1063 | != travelling from the top of the atmosphere to layer i; =* |
---|
1064 | != NID(i), i = 0..NZ-1 =* |
---|
1065 | != DELTA - LOGICAL, switch to use delta-scaling (I)=* |
---|
1066 | != .TRUE. -> apply delta-scaling =* |
---|
1067 | != .FALSE.-> do not apply delta-scaling =* |
---|
1068 | != FDR - REAL, contribution of the direct component to the total (O)=* |
---|
1069 | != actinic flux at each altitude level =* |
---|
1070 | != FUP - REAL, contribution of the diffuse upwelling component to (O)=* |
---|
1071 | != the total actinic flux at each altitude level =* |
---|
1072 | != FDN - REAL, contribution of the diffuse downwelling component to (O)=* |
---|
1073 | != the total actinic flux at each altitude level =* |
---|
1074 | != EDR - REAL, contribution of the direct component to the total (O)=* |
---|
1075 | != spectral irradiance at each altitude level =* |
---|
1076 | != EUP - REAL, contribution of the diffuse upwelling component to (O)=* |
---|
1077 | != the total spectral irradiance at each altitude level =* |
---|
1078 | != EDN - REAL, contribution of the diffuse downwelling component to (O)=* |
---|
1079 | *= the total spectral irradiance at each altitude level =* |
---|
1080 | !-----------------------------------------------------------------------------* |
---|
1081 | |
---|
1082 | implicit none |
---|
1083 | |
---|
1084 | ! input: |
---|
1085 | |
---|
1086 | INTEGER nlev |
---|
1087 | REAL zen, rsfc |
---|
1088 | REAL tauu(nlev), omu(nlev), gu(nlev) |
---|
1089 | REAL dsdh(0:nlev,nlev) |
---|
1090 | INTEGER nid(0:nlev) |
---|
1091 | LOGICAL delta |
---|
1092 | |
---|
1093 | ! output: |
---|
1094 | |
---|
1095 | REAL fup(nlev),fdn(nlev),fdr(nlev) |
---|
1096 | REAL eup(nlev),edn(nlev),edr(nlev) |
---|
1097 | |
---|
1098 | ! local: |
---|
1099 | |
---|
1100 | REAL tausla(0:nlev), tauc(0:nlev) |
---|
1101 | REAL mu2(0:nlev), mu, sum |
---|
1102 | |
---|
1103 | ! internal coefficients and matrix |
---|
1104 | |
---|
1105 | REAL lam(nlev),taun(nlev),bgam(nlev) |
---|
1106 | REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev) |
---|
1107 | REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev) |
---|
1108 | REAL mu1(nlev) |
---|
1109 | INTEGER row |
---|
1110 | REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev) |
---|
1111 | |
---|
1112 | ! other: |
---|
1113 | |
---|
1114 | REAL pifs, fdn0 |
---|
1115 | REAL gi(nlev), omi(nlev), tempg |
---|
1116 | REAL f, g, om |
---|
1117 | REAL gam1, gam2, gam3, gam4 |
---|
1118 | real, parameter :: largest = 1.e+36 |
---|
1119 | real, parameter :: precis = 1.e-7 |
---|
1120 | |
---|
1121 | ! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4 |
---|
1122 | ! in delta-function, modified quadrature, hemispheric constant, |
---|
1123 | ! Hybrid modified Eddington-delta function metods, p633,Table1. |
---|
1124 | ! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 |
---|
1125 | ! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, |
---|
1126 | ! uncomment the following two lines and the appropriate statements further |
---|
1127 | ! down. |
---|
1128 | ! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0, |
---|
1129 | ! > BETA1, BETAn, amu1, subd |
---|
1130 | |
---|
1131 | REAL expon, expon0, expon1, divisr, temp, up, dn |
---|
1132 | REAL ssfc |
---|
1133 | INTEGER nlayer, mrows, lev |
---|
1134 | |
---|
1135 | INTEGER i, j |
---|
1136 | |
---|
1137 | ! Some additional program constants: |
---|
1138 | |
---|
1139 | real pi, dr |
---|
1140 | REAL eps |
---|
1141 | PARAMETER (eps = 1.E-3) |
---|
1142 | !_______________________________________________________________________ |
---|
1143 | |
---|
1144 | ! MU = cosine of solar zenith angle |
---|
1145 | ! RSFC = surface albedo |
---|
1146 | ! TAUU = unscaled optical depth of each layer |
---|
1147 | ! OMU = unscaled single scattering albedo |
---|
1148 | ! GU = unscaled asymmetry factor |
---|
1149 | ! KLEV = max dimension of number of layers in atmosphere |
---|
1150 | ! NLAYER = number of layers in the atmosphere |
---|
1151 | ! NLEVEL = nlayer + 1 = number of levels |
---|
1152 | |
---|
1153 | ! initial conditions: pi*solar flux = 1; diffuse incidence = 0 |
---|
1154 | |
---|
1155 | pifs = 1. |
---|
1156 | fdn0 = 0. |
---|
1157 | |
---|
1158 | nlayer = nlev - 1 |
---|
1159 | |
---|
1160 | pi = acos(-1.) |
---|
1161 | dr = pi/180. |
---|
1162 | mu = COS(zen*dr) |
---|
1163 | |
---|
1164 | !************* compute coefficients for each layer: |
---|
1165 | ! GAM1 - GAM4 = 2-stream coefficients, different for different approximations |
---|
1166 | ! EXPON0 = calculation of e when TAU is zero |
---|
1167 | ! EXPON1 = calculation of e when TAU is TAUN |
---|
1168 | ! CUP and CDN = calculation when TAU is zero |
---|
1169 | ! CUPTN and CDNTN = calc. when TAU is TAUN |
---|
1170 | ! DIVISR = prevents division by zero |
---|
1171 | |
---|
1172 | do j = 0, nlev |
---|
1173 | tauc(j) = 0. |
---|
1174 | tausla(j) = 0. |
---|
1175 | mu2(j) = 1./SQRT(largest) |
---|
1176 | end do |
---|
1177 | |
---|
1178 | IF (.NOT. delta) THEN |
---|
1179 | DO i = 1, nlayer |
---|
1180 | gi(i) = gu(i) |
---|
1181 | omi(i) = omu(i) |
---|
1182 | taun(i) = tauu(i) |
---|
1183 | END DO |
---|
1184 | ELSE |
---|
1185 | |
---|
1186 | ! delta-scaling. Have to be done for delta-Eddington approximation, |
---|
1187 | ! delta discrete ordinate, Practical Improved Flux Method, delta function, |
---|
1188 | ! and Hybrid modified Eddington-delta function methods approximations |
---|
1189 | |
---|
1190 | DO i = 1, nlayer |
---|
1191 | f = gu(i)*gu(i) |
---|
1192 | gi(i) = (gu(i) - f)/(1 - f) |
---|
1193 | omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f) |
---|
1194 | taun(i) = (1 - omu(i)*f)*tauu(i) |
---|
1195 | END DO |
---|
1196 | END IF |
---|
1197 | |
---|
1198 | ! calculate slant optical depth at the top of the atmosphere when zen>90. |
---|
1199 | ! in this case, higher altitude of the top layer is recommended. |
---|
1200 | |
---|
1201 | IF (zen .GT. 90.0) THEN |
---|
1202 | IF (nid(0) .LT. 0) THEN |
---|
1203 | tausla(0) = largest |
---|
1204 | ELSE |
---|
1205 | sum = 0.0 |
---|
1206 | DO j = 1, nid(0) |
---|
1207 | sum = sum + 2.*taun(j)*dsdh(0,j) |
---|
1208 | END DO |
---|
1209 | tausla(0) = sum |
---|
1210 | END IF |
---|
1211 | END IF |
---|
1212 | |
---|
1213 | DO 11, i = 1, nlayer |
---|
1214 | g = gi(i) |
---|
1215 | om = omi(i) |
---|
1216 | tauc(i) = tauc(i-1) + taun(i) |
---|
1217 | |
---|
1218 | ! stay away from 1 by precision. For g, also stay away from -1 |
---|
1219 | |
---|
1220 | tempg = AMIN1(abs(g),1. - precis) |
---|
1221 | g = SIGN(tempg,g) |
---|
1222 | om = AMIN1(om,1.-precis) |
---|
1223 | |
---|
1224 | ! calculate slant optical depth |
---|
1225 | |
---|
1226 | IF (nid(i) .LT. 0) THEN |
---|
1227 | tausla(i) = largest |
---|
1228 | ELSE |
---|
1229 | sum = 0.0 |
---|
1230 | DO j = 1, MIN(nid(i),i) |
---|
1231 | sum = sum + taun(j)*dsdh(i,j) |
---|
1232 | END DO |
---|
1233 | DO j = MIN(nid(i),i)+1,nid(i) |
---|
1234 | sum = sum + 2.*taun(j)*dsdh(i,j) |
---|
1235 | END DO |
---|
1236 | tausla(i) = sum |
---|
1237 | IF (tausla(i) .EQ. tausla(i-1)) THEN |
---|
1238 | mu2(i) = SQRT(largest) |
---|
1239 | ELSE |
---|
1240 | mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1)) |
---|
1241 | mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)), |
---|
1242 | $ mu2(i) ) |
---|
1243 | END IF |
---|
1244 | END IF |
---|
1245 | |
---|
1246 | !** the following gamma equations are from pg 16,289, Table 1 |
---|
1247 | !** save mu1 for each approx. for use in converting irradiance to actinic flux |
---|
1248 | |
---|
1249 | ! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452): |
---|
1250 | |
---|
1251 | c gam1 = (7. - om*(4. + 3.*g))/4. |
---|
1252 | c gam2 = -(1. - om*(4. - 3.*g))/4. |
---|
1253 | c gam3 = (2. - 3.*g*mu)/4. |
---|
1254 | c gam4 = 1. - gam3 |
---|
1255 | c mu1(i) = 0.5 |
---|
1256 | |
---|
1257 | * quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475): |
---|
1258 | |
---|
1259 | c gam1 = 1.7320508*(2. - om*(1. + g))/2. |
---|
1260 | c gam2 = 1.7320508*om*(1. - g)/2. |
---|
1261 | c gam3 = (1. - 1.7320508*g*mu)/2. |
---|
1262 | c gam4 = 1. - gam3 |
---|
1263 | c mu1(i) = 1./sqrt(3.) |
---|
1264 | |
---|
1265 | * hemispheric mean (Toon et al., 1089, JGR, 94, 16287): |
---|
1266 | |
---|
1267 | gam1 = 2. - om*(1. + g) |
---|
1268 | gam2 = om*(1. - g) |
---|
1269 | gam3 = (2. - g*mu)/4. |
---|
1270 | gam4 = 1. - gam3 |
---|
1271 | mu1(i) = 0.5 |
---|
1272 | |
---|
1273 | * PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166): |
---|
1274 | c GAM1 = 0.25*(8. - OM*(5. + 3.*G)) |
---|
1275 | c GAM2 = 0.75*OM*(1.-G) |
---|
1276 | c GAM3 = 0.25*(2.-3.*G*MU) |
---|
1277 | c GAM4 = 1. - GAM3 |
---|
1278 | c mu1(i) = 0.5 |
---|
1279 | |
---|
1280 | * delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26): |
---|
1281 | c GAM1 = 0.5*1.7320508*(2. - OM*(1. + G)) |
---|
1282 | c GAM2 = 0.5*1.7320508*OM*(1.-G) |
---|
1283 | c GAM3 = 0.5*(1.-1.7320508*G*MU) |
---|
1284 | c GAM4 = 1. - GAM3 |
---|
1285 | c mu1(i) = 1./sqrt(3.) |
---|
1286 | |
---|
1287 | * Calculations of Associated Legendre Polynomials for GAMA1,2,3,4 |
---|
1288 | * in delta-function, modified quadrature, hemispheric constant, |
---|
1289 | * Hybrid modified Eddington-delta function metods, p633,Table1. |
---|
1290 | * W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630 |
---|
1291 | * W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440 |
---|
1292 | c YLM0 = 2. |
---|
1293 | c YLM2 = -3.*G*MU |
---|
1294 | c YLM4 = 0.875*G**3*MU*(5.*MU**2-3.) |
---|
1295 | c YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4) |
---|
1296 | c YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4 |
---|
1297 | c *+429.*MU**6) |
---|
1298 | c YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4 |
---|
1299 | c *-25740.*MU**6+12155.*MU**8) |
---|
1300 | c YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4 |
---|
1301 | c *+218790.*MU**6-230945.*MU**8+88179.*MU**10) |
---|
1302 | c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12 |
---|
1303 | c YLMS=0.25*YLMS |
---|
1304 | c BETA0 = YLMS |
---|
1305 | c |
---|
1306 | c amu1=1./1.7320508 |
---|
1307 | c YLM0 = 2. |
---|
1308 | c YLM2 = -3.*G*amu1 |
---|
1309 | c YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.) |
---|
1310 | c YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4) |
---|
1311 | c YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4 |
---|
1312 | c *+429.*amu1**6) |
---|
1313 | c YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4 |
---|
1314 | c *-25740.*amu1**6+12155.*amu1**8) |
---|
1315 | c YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4 |
---|
1316 | c *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10) |
---|
1317 | c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12 |
---|
1318 | c YLMS=0.25*YLMS |
---|
1319 | c BETA1 = YLMS |
---|
1320 | c |
---|
1321 | c BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5 |
---|
1322 | c *-0.045776*G**7) |
---|
1323 | |
---|
1324 | |
---|
1325 | * Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630): |
---|
1326 | c subd=4.*(1.-G*G*(1.-MU)) |
---|
1327 | c GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd |
---|
1328 | c GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd |
---|
1329 | c GAM3 = BETA0 |
---|
1330 | c GAM4 = 1. - GAM3 |
---|
1331 | c mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g) |
---|
1332 | |
---|
1333 | ***** |
---|
1334 | * delta function (Meador, and Weaver, 1980, JAS, 37, 630): |
---|
1335 | c GAM1 = (1. - OM*(1. - beta0))/MU |
---|
1336 | c GAM2 = OM*BETA0/MU |
---|
1337 | c GAM3 = BETA0 |
---|
1338 | c GAM4 = 1. - GAM3 |
---|
1339 | c mu1(i) = mu |
---|
1340 | ***** |
---|
1341 | * modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630): |
---|
1342 | c GAM1 = 1.7320508*(1. - OM*(1. - beta1)) |
---|
1343 | c GAM2 = 1.7320508*OM*beta1 |
---|
1344 | c GAM3 = BETA0 |
---|
1345 | c GAM4 = 1. - GAM3 |
---|
1346 | c mu1(i) = 1./sqrt(3.) |
---|
1347 | |
---|
1348 | * hemispheric constant (Toon et al., 1989, JGR, 94, 16287): |
---|
1349 | c GAM1 = 2.*(1. - OM*(1. - betan)) |
---|
1350 | c GAM2 = 2.*OM*BETAn |
---|
1351 | c GAM3 = BETA0 |
---|
1352 | c GAM4 = 1. - GAM3 |
---|
1353 | c mu1(i) = 0.5 |
---|
1354 | |
---|
1355 | ***** |
---|
1356 | |
---|
1357 | * lambda = pg 16,290 equation 21 |
---|
1358 | * big gamma = pg 16,290 equation 22 |
---|
1359 | * if gam2 = 0., then bgam = 0. |
---|
1360 | |
---|
1361 | lam(i) = sqrt(gam1*gam1 - gam2*gam2) |
---|
1362 | |
---|
1363 | IF (gam2 .NE. 0.) THEN |
---|
1364 | bgam(i) = (gam1 - lam(i))/gam2 |
---|
1365 | ELSE |
---|
1366 | bgam(i) = 0. |
---|
1367 | END IF |
---|
1368 | |
---|
1369 | expon = EXP(-lam(i)*taun(i)) |
---|
1370 | |
---|
1371 | * e1 - e4 = pg 16,292 equation 44 |
---|
1372 | |
---|
1373 | e1(i) = 1. + bgam(i)*expon |
---|
1374 | e2(i) = 1. - bgam(i)*expon |
---|
1375 | e3(i) = bgam(i) + expon |
---|
1376 | e4(i) = bgam(i) - expon |
---|
1377 | |
---|
1378 | * the following sets up for the C equations 23, and 24 |
---|
1379 | * found on page 16,290 |
---|
1380 | * prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3 |
---|
1381 | * which is approx equiv to shifting MU by 0.5*EPS* (MU)**3 |
---|
1382 | |
---|
1383 | expon0 = EXP(-tausla(i-1)) |
---|
1384 | expon1 = EXP(-tausla(i)) |
---|
1385 | |
---|
1386 | divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i)) |
---|
1387 | temp = AMAX1(eps,abs(divisr)) |
---|
1388 | divisr = SIGN(temp,divisr) |
---|
1389 | |
---|
1390 | up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr |
---|
1391 | dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr |
---|
1392 | |
---|
1393 | * cup and cdn are when tau is equal to zero |
---|
1394 | * cuptn and cdntn are when tau is equal to taun |
---|
1395 | |
---|
1396 | cup(i) = up*expon0 |
---|
1397 | cdn(i) = dn*expon0 |
---|
1398 | cuptn(i) = up*expon1 |
---|
1399 | cdntn(i) = dn*expon1 |
---|
1400 | |
---|
1401 | 11 CONTINUE |
---|
1402 | |
---|
1403 | ***************** set up matrix ****** |
---|
1404 | * ssfc = pg 16,292 equation 37 where pi Fs is one (unity). |
---|
1405 | |
---|
1406 | ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs |
---|
1407 | |
---|
1408 | * MROWS = the number of rows in the matrix |
---|
1409 | |
---|
1410 | mrows = 2*nlayer |
---|
1411 | |
---|
1412 | * the following are from pg 16,292 equations 39 - 43. |
---|
1413 | * set up first row of matrix: |
---|
1414 | |
---|
1415 | i = 1 |
---|
1416 | a(1) = 0. |
---|
1417 | b(1) = e1(i) |
---|
1418 | d(1) = -e2(i) |
---|
1419 | e(1) = fdn0 - cdn(i) |
---|
1420 | |
---|
1421 | row=1 |
---|
1422 | |
---|
1423 | * set up odd rows 3 thru (MROWS - 1): |
---|
1424 | |
---|
1425 | i = 0 |
---|
1426 | DO 20, row = 3, mrows - 1, 2 |
---|
1427 | i = i + 1 |
---|
1428 | a(row) = e2(i)*e3(i) - e4(i)*e1(i) |
---|
1429 | b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1) |
---|
1430 | d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1) |
---|
1431 | e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + |
---|
1432 | $ e1(i)*(cdntn(i) - cdn(i + 1)) |
---|
1433 | 20 CONTINUE |
---|
1434 | |
---|
1435 | * set up even rows 2 thru (MROWS - 2): |
---|
1436 | |
---|
1437 | i = 0 |
---|
1438 | DO 30, row = 2, mrows - 2, 2 |
---|
1439 | i = i + 1 |
---|
1440 | a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1) |
---|
1441 | b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1) |
---|
1442 | d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1) |
---|
1443 | e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - |
---|
1444 | $ (cdn(i + 1) - cdntn(i))*e4(i + 1) |
---|
1445 | 30 CONTINUE |
---|
1446 | |
---|
1447 | * set up last row of matrix at MROWS: |
---|
1448 | |
---|
1449 | row = mrows |
---|
1450 | i = nlayer |
---|
1451 | |
---|
1452 | a(row) = e1(i) - rsfc*e3(i) |
---|
1453 | b(row) = e2(i) - rsfc*e4(i) |
---|
1454 | d(row) = 0. |
---|
1455 | e(row) = ssfc - cuptn(i) + rsfc*cdntn(i) |
---|
1456 | |
---|
1457 | * solve tri-diagonal matrix: |
---|
1458 | |
---|
1459 | CALL tridiag(a, b, d, e, y, mrows) |
---|
1460 | |
---|
1461 | **** unfold solution of matrix, compute output fluxes: |
---|
1462 | |
---|
1463 | row = 1 |
---|
1464 | lev = 1 |
---|
1465 | j = 1 |
---|
1466 | |
---|
1467 | * the following equations are from pg 16,291 equations 31 & 32 |
---|
1468 | |
---|
1469 | fdr(lev) = EXP( -tausla(0) ) |
---|
1470 | edr(lev) = mu * fdr(lev) |
---|
1471 | edn(lev) = fdn0 |
---|
1472 | eup(lev) = y(row)*e3(j) - y(row + 1)*e4(j) + cup(j) |
---|
1473 | fdn(lev) = edn(lev)/mu1(lev) |
---|
1474 | fup(lev) = eup(lev)/mu1(lev) |
---|
1475 | |
---|
1476 | DO 60, lev = 2, nlayer + 1 |
---|
1477 | fdr(lev) = EXP(-tausla(lev-1)) |
---|
1478 | edr(lev) = mu *fdr(lev) |
---|
1479 | edn(lev) = y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j) |
---|
1480 | eup(lev) = y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j) |
---|
1481 | fdn(lev) = edn(lev)/mu1(j) |
---|
1482 | fup(lev) = eup(lev)/mu1(j) |
---|
1483 | |
---|
1484 | row = row + 2 |
---|
1485 | j = j + 1 |
---|
1486 | 60 CONTINUE |
---|
1487 | |
---|
1488 | end subroutine ps2str |
---|
1489 | |
---|
1490 | *=============================================================================* |
---|
1491 | |
---|
1492 | subroutine tridiag(a,b,c,r,u,n) |
---|
1493 | |
---|
1494 | !_______________________________________________________________________ |
---|
1495 | ! solves tridiagonal system. From Numerical Recipies, p. 40 |
---|
1496 | !_______________________________________________________________________ |
---|
1497 | |
---|
1498 | IMPLICIT NONE |
---|
1499 | |
---|
1500 | ! input: |
---|
1501 | |
---|
1502 | INTEGER n |
---|
1503 | REAL a, b, c, r |
---|
1504 | DIMENSION a(n),b(n),c(n),r(n) |
---|
1505 | |
---|
1506 | ! output: |
---|
1507 | |
---|
1508 | REAL u |
---|
1509 | DIMENSION u(n) |
---|
1510 | |
---|
1511 | ! local: |
---|
1512 | |
---|
1513 | INTEGER j |
---|
1514 | |
---|
1515 | REAL bet, gam |
---|
1516 | DIMENSION gam(n) |
---|
1517 | !_______________________________________________________________________ |
---|
1518 | |
---|
1519 | IF (b(1) .EQ. 0.) STOP 1001 |
---|
1520 | bet = b(1) |
---|
1521 | u(1) = r(1)/bet |
---|
1522 | DO 11, j = 2, n |
---|
1523 | gam(j) = c(j - 1)/bet |
---|
1524 | bet = b(j) - a(j)*gam(j) |
---|
1525 | IF (bet .EQ. 0.) STOP 2002 |
---|
1526 | u(j) = (r(j) - a(j)*u(j - 1))/bet |
---|
1527 | 11 CONTINUE |
---|
1528 | DO 12, j = n - 1, 1, -1 |
---|
1529 | u(j) = u(j) - gam(j + 1)*u(j + 1) |
---|
1530 | 12 CONTINUE |
---|
1531 | !_______________________________________________________________________ |
---|
1532 | |
---|
1533 | end subroutine tridiag |
---|
1534 | |
---|
1535 | end subroutine photolysis_online |
---|