Last change
on this file since 943 was
943,
checked in by aslmd, 12 years ago
|
clean and organized UTIL/PYTHON folder
|
-
Property svn:executable set to
*
|
File size:
1.6 KB
|
Line | |
---|
1 | #! /usr/bin/env python |
---|
2 | |
---|
3 | import numpy as np |
---|
4 | import matplotlib.pyplot as mpl |
---|
5 | from scipy import integrate |
---|
6 | |
---|
7 | def singlet(rel_rho,z): |
---|
8 | |
---|
9 | ###### z MUST BE IN km |
---|
10 | z = np.array(z) |
---|
11 | rel_rho = np.array(rel_rho) |
---|
12 | |
---|
13 | alpha = 1.7e-3 ##sec^(-1)= ozone photolysis rate constant - from Nair et al. 1994, scaled at the considered LS |
---|
14 | k = 2.0e-20 ##cm^(-3)molecules^(-1) sec^(-1) = co2 quenching rate constant - upper limit from |
---|
15 | tau = 3880 ##sec = O2(a1Delatg) molecules life time |
---|
16 | gamma = 1.4 ##;= ratio of specific heats |
---|
17 | zh_o3 = 10 ##; *** km = peak altitude of the O3 layer |
---|
18 | H = 11.0 ##;scale eight |
---|
19 | H = 8.3 |
---|
20 | H = 9.5 |
---|
21 | co2_0 = 2e17 |
---|
22 | exponent = - z / H |
---|
23 | co2p_u = co2_0 * np.exp(exponent) |
---|
24 | #;analytical expression for the O3 unperturbed vertical profile |
---|
25 | zo3 = 26. |
---|
26 | H1 = 5.4 |
---|
27 | H2 = 15. |
---|
28 | f = 1.17e9 |
---|
29 | |
---|
30 | #### COMPUTE g0 |
---|
31 | exponent = -(z-zo3)/H2 |
---|
32 | g0 = (gamma*H-H1)/((gamma-1.)*H1) -(gamma*H)/(H2*(gamma-1.))*np.exp(exponent) |
---|
33 | |
---|
34 | ### OZONE PERTURBED AND UNPERTURBED |
---|
35 | exponent1 = -(z-zo3)/H1 |
---|
36 | exponent2 = -(z-zo3)/H2 |
---|
37 | o3p_u = f * np.exp( exponent1 - np.exp(exponent2) ) |
---|
38 | o3p_p = o3p_u*np.power(rel_rho,g0) |
---|
39 | |
---|
40 | ### CO2 STUFF |
---|
41 | co2p_p = co2p_u*rel_rho |
---|
42 | |
---|
43 | ### SINGLET PROFILE |
---|
44 | o2p_p = (0.9*alpha*o3p_p)/(1.+tau*k*co2p_p) |
---|
45 | o2p_u = (0.9*alpha*o3p_u)/(1.+tau*k*co2p_u) |
---|
46 | rel = 100. * ( o2p_p / o2p_u - 1. ) |
---|
47 | |
---|
48 | ### SINGLET INTEGRATED |
---|
49 | o2p_pi = integrate.trapz(o2p_p,x=z) |
---|
50 | o2p_ui = integrate.trapz(o2p_u,x=z) |
---|
51 | reli = 100. * ( o2p_pi / o2p_ui - 1. ) |
---|
52 | |
---|
53 | #mpl.plot(rel,z) |
---|
54 | #mpl.show() |
---|
55 | |
---|
56 | return reli, rel |
---|
Note: See
TracBrowser
for help on using the repository browser.