source: trunk/UTIL/PYTHON/singlet.py @ 805

#! /usr/bin/env python

import numpy as np
import matplotlib.pyplot as mpl
from scipy import integrate 

def singlet(rel_rho,z):

        ###### z MUST BE IN km
        z = np.array(z)
        rel_rho = np.array(rel_rho)

        alpha = 1.7e-3          ##sec^(-1)= ozone photolysis rate constant - from Nair et al. 1994, scaled at the considered LS
        k =  2.0e-20            ##cm^(-3)molecules^(-1) sec^(-1) = co2 quenching rate constant - upper limit from 
        tau = 3880              ##sec = O2(a1Delatg) molecules life time
        gamma = 1.4             ##;= ratio of specific heats
        zh_o3 = 10              ##; *** km = peak altitude of the O3 layer 
        H = 11.0                ##;scale eight
        H = 8.3
        H = 9.5
        co2_0 = 2e17
        exponent = - z / H
        co2p_u = co2_0 * np.exp(exponent) 
        #;analytical expression for the O3 unperturbed vertical profile 
        zo3 = 26.
        H1 = 5.4
        H2 = 15.
        f = 1.17e9

        #### COMPUTE g0
        exponent = -(z-zo3)/H2
        g0 = (gamma*H-H1)/((gamma-1.)*H1) -(gamma*H)/(H2*(gamma-1.))*np.exp(exponent)

        ### OZONE PERTURBED AND UNPERTURBED
        exponent1 = -(z-zo3)/H1
        exponent2 = -(z-zo3)/H2
        o3p_u = f * np.exp( exponent1 - np.exp(exponent2) )
        o3p_p = o3p_u*np.power(rel_rho,g0)

        ### CO2 STUFF
        co2p_p = co2p_u*rel_rho

        ### SINGLET PROFILE
        o2p_p = (0.9*alpha*o3p_p)/(1.+tau*k*co2p_p)     
        o2p_u = (0.9*alpha*o3p_u)/(1.+tau*k*co2p_u) 
        rel  = 100. * ( o2p_p  / o2p_u  - 1. )

        ### SINGLET INTEGRATED
        o2p_pi = integrate.trapz(o2p_p,x=z)
        o2p_ui = integrate.trapz(o2p_u,x=z)
        reli = 100. * ( o2p_pi / o2p_ui - 1. )

        #mpl.plot(rel,z)
        #mpl.show()

        return reli, rel