#!/usr/bin/env python
# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt

# LIQUID PRECIP ---------------------------------------------------

qlw = 0.1e-3 # kg/kg
tau = 900. # seconds
clw = 0.65e-3 # kg/kg

dqlw = -qlw/tau*(1.-np.exp(-(qlw/clw)**2.))

var=dqlw*3600.*1e3 # g/kg/hr
print("dqlw=",var,"g/kg/hr")

# SOLID PRECIP ----------------------------------------------------

rho=750e2/287./273.15 # kg/m3
wiw=10e-2 # m/s
qiw_zB = 0.0e-3 # kg/kg
qiw_zA = 0.1e-3 # kg/kg
zB = 3000. # m
zA = 2500. # m
mflux_zA = rho * wiw * qiw_zA
mflux_zB = rho * wiw * qiw_zB

dqiw = 1./rho * (mflux_zB-mflux_zA) / (zB-zA)

var=dqiw*3600.*1e3 # g/kg/hour
print("dqiw=",var,"g/kg/hr")


# AUTOCONVERSION --------------------------------------------------

fig = plt.figure()

tau = 1800. # seconds
qlw = np.linspace(0,1e-3)
clw = 0.65e-3 # kg/kg
dqlw = -qlw/tau*(1.-np.exp(-(qlw/clw)**2.))
plt.plot(qlw*1e3, dqlw*3600.*1e3, 'b-',linewidth=2,
        label=r"$\tau$=1800s, clw=0.65 g/kg")
plt.gca().set_xlabel(r'q$_{l} (g/kg)$')
plt.gca().set_ylabel(r'$d_{t}q_{l} (g/kg/hr)$')

tau = 900. # seconds
dqlw = -qlw/tau*(1.-np.exp(-(qlw/clw)**2.))
plt.plot(qlw*1e3, dqlw*3600.*1e3, 'r-',linewidth=2,
        label=r"$\tau$=900s, clw=0.65 g/kg")

tau = 600. # seconds
dqlw = -qlw/tau*(1.-np.exp(-(qlw/clw)**2.))
plt.plot(qlw*1e3, dqlw*3600.*1e3, 'g-',linewidth=2,
        label=r"$\tau$=600s, clw=0.65 g/kg")

tau = 600. # seconds
clw = 0.4e-3 # kg/kg
dqlw = -qlw/tau*(1.-np.exp(-(qlw/clw)**2.))
plt.plot(qlw*1e3, dqlw*3600.*1e3, 'g--',linewidth=2,
        label=r"$\tau$=600s, clw=0.4 g/kg")

plt.legend(fontsize=11,loc='best')

plt.show()