Main model
This “main” program imported all the modules, and ran the total time for-loop.
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.basemap import Basemap
from parameters import *
from variables import *
from initial import *
from advectionx import*
from advectiony import*
from advectionz import*
from dispersionx import *
from dispersiony import *
from dispersionz import *
from degree_conversion import *
from settling import *
from landed import *
#[height,lat,lon]; [z,y,x]
Rold = initial(sigma_lat,sigma_lon)
lon1 = 65.
lon2 = 155.
lat1 = -25.
lat2 = 25.
fig = plt.figure(1,facecolor = 'white')
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('longitude [decimal degees]')
ax.set_ylabel('latitude [decimal degees]')
ax.set_zlabel('altitude [meters]')
ax.set_zlim(0,20000)
ax.set_ylim(lat1,lat2)
ax.set_xlim(lon1,lon2)
ax.scatter(Rold[:,2],Rold[:,1],Rold[:,0],c='r',s=10)
ax.set_title('Particle Trajectory')
plt.figure(2)
m = Basemap(llcrnrlon=lon1,llcrnrlat=lat1,urcrnrlon=lon2,urcrnrlat=lat2,
projection='lcc',lon_0=start_lon,lat_0=start_lat,
resolution ='l',area_thresh=1000.)
m.drawcoastlines(linewidth=0.25)
m.drawcountries(linewidth=0.25)
m.drawparallels(np.arange(lat1,lat2,10),labels=[1,1,0,0])
m.drawmeridians(np.arange(lon1,lon2,10),labels=[0,0,0,1])
x, y = m(Rold[:,2],Rold[:,1])
m.scatter(x,y,s=20,marker='o',color='r',zorder=10)
plt.title('Particle Trajectory')
j=3
for t in xrange(1,time_steps):
t_index = int(float(t)/float(time_steps)*20.)
t_time_1 = float(t)/float(time_steps)*20.
t_time_2 = (float(t)+0.5)/float(time_steps)*20.
t_time_3 = (float(t)+0.5)/float(time_steps)*20.
t_time_4 = (float(t)+1.0)/float(time_steps)*20.
dx_lat, dx_lon = degree_to_dx(R)
#RUNGA KUTTA
R_step1[:,:] = 0.0
R_step2[:,:] = 0.0
R_step3[:,:] = 0.0
R_step4[:,:] = 0.0
#dispersion
dispz_temp = dispersionz(Active) * dt
dispy_temp = dispersiony(Active) * dt / dx_lat
dispx_temp = dispersionx(Active) * dt / dx_lon
R_step1[:,0] += dispz_temp
R_step2[:,0] += dispz_temp
R_step3[:,0] += dispz_temp
R_step4[:,0] += dispz_temp
R_step1[:,1] += dispy_temp
R_step2[:,1] += dispy_temp
R_step3[:,1] += dispy_temp
R_step4[:,1] += dispy_temp
R_step1[:,2] += dispx_temp
R_step2[:,2] += dispx_temp
R_step3[:,2] += dispx_temp
R_step4[:,2] += dispx_temp
#settling
R_step1[:,0] -= settling(D,Active) * dt
R_step2[:,0] -= settling(D,Active) * dt
R_step3[:,0] -= settling(D,Active) * dt
R_step4[:,0] -= settling(D,Active) * dt
R_step1[:,0] += advectionz(Rold,Active,t_index,t_time_1) * dt
R_step1[:,1] += advectiony(Rold,Active,t_index,t_time_1) * dt / dx_lat
R_step1[:,2] += advectionx(Rold,Active,t_index,t_time_1) * dt / dx_lon
R_step2[:,0] += advectionz(Rold+R_step1/2.0,Active,t_index,t_time_2) * dt
R_step2[:,1] += advectiony(Rold+R_step1/2.0,Active,t_index,t_time_2) * dt / dx_lat
R_step2[:,2] += advectionx(Rold+R_step1/2.0,Active,t_index,t_time_2) * dt / dx_lon
R_step3[:,0] += advectionz(Rold+R_step2/2.0,Active,t_index,t_time_3) * dt
R_step3[:,1] += advectiony(Rold+R_step2/2.0,Active,t_index,t_time_3) * dt / dx_lat
R_step3[:,2] += advectionx(Rold+R_step2/2.0,Active,t_index,t_time_3) * dt / dx_lon
R_step4[:,0] += advectionz(Rold+R_step3,Active,t_index,t_time_4) * dt
R_step4[:,1] += advectiony(Rold+R_step3,Active,t_index,t_time_4) * dt / dx_lat
R_step4[:,2] += advectionx(Rold+R_step3,Active,t_index,t_time_4) * dt / dx_lon
Rnew = Rold + R_step1/6.0 + R_step2/3.0 + R_step3/3.0 + R_step4/6.0
Active, Rnew[:,0] = landed(Rnew[:,0],Active)
if t%(time_steps/50)==0:
ax.scatter(Rnew[:,2],Rnew[:,1],Rnew[:,0],c='k',s=1,alpha=float(t)/float(time_steps-1)*0.5)
plt.figure(2,facecolor='white')
x, y = m(Rnew[:,2],Rnew[:,1])
m.scatter(x,y,s=1,marker='o',color='k',alpha=float(t)/float(time_steps-1)*0.5,zorder=10)
print str(int(100.*float(t)/float(time_steps-1))) + '% done'
if t%(time_steps/10)==0:
plt.figure(j)
m = Basemap(llcrnrlon=lon1,llcrnrlat=lat1,urcrnrlon=lon2,urcrnrlat=lat2,
projection='lcc',lon_0=start_lon,lat_0=start_lat,
resolution ='l',area_thresh=1000.)
m.drawcoastlines(linewidth=0.25)
m.drawcountries(linewidth=0.25)
m.drawparallels(np.arange(lat1,lat2,10),labels=[1,1,0,0])
m.drawmeridians(np.arange(lon1,lon2,10),labels=[0,0,0,1])
x, y = m(start_lon,start_lat)
m.scatter(x,y,s=20,marker='o',color='r',zorder=10)
x, y = m(Rnew[:,2],Rnew[:,1])
m.scatter(x,y,s=2,marker='o',color='k',zorder=10,alpha=0.5)
plt.title('Plume after ' + str((j-2)*12) + ' hrs')
fig = plt.figure(j+10,facecolor = 'white')
axx = fig.add_subplot(111, projection='3d')
axx.set_xlabel('longitude [decimal degees]')
axx.set_ylabel('latitude [decimal degees]')
axx.set_zlabel('altitude [meters]')
axx.set_zlim(0,20000)
axx.set_ylim(lat1,lat2)
axx.set_xlim(lon1,lon2)
axx.scatter(Rnew[:,2],Rnew[:,1],Rnew[:,0],c='k',s=1)
axx.set_title('Plume after ' + str((j-2)*12) + ' hrs')
j+=1
Rold = Rnew
for i in xrange(0,num_particles):
if Active[i] == 0:
ax.scatter(Rnew[i,2],Rnew[i,1],Rnew[i,0],c='g',s=10)
plt.figure(2,facecolor='white')
x, y = m(Rnew[i,2],Rnew[i,1])
m.scatter(x,y,s=2,marker='o',color='g',zorder=10)
else:
ax.scatter(Rnew[i,2],Rnew[i,1],Rnew[i,0],c='b',s=10)
plt.figure(2,facecolor='white')
x, y = m(Rnew[i,2],Rnew[i,1])
m.scatter(x,y,s=10,marker='o',color='b',zorder=10)
plt.figure(1,facecolor='white')
plt.savefig('3D-final.png')
plt.figure(2,facecolor='white')
plt.savefig('2D-final.png')
plt.figure(3,facecolor='white')
plt.savefig('2D-12hr.png')
plt.figure(4,facecolor='white')
plt.savefig('2D-24hr.png')
plt.figure(5,facecolor='white')
plt.savefig('2D-36hr.png')
plt.figure(6,facecolor='white')
plt.savefig('2D-48hr.png')
plt.figure(7,facecolor='white')
plt.savefig('2D-60hr.png')
plt.figure(8,facecolor='white')
plt.savefig('2D-72hr.png')
plt.figure(9,facecolor='white')
plt.savefig('2D-84hr.png')
plt.figure(10,facecolor='white')
plt.savefig('2D-96hr.png')
plt.figure(11,facecolor='white')
plt.savefig('2D-108hr.png')
plt.figure(12,facecolor='white')
plt.savefig('2D-120hr.png')
plt.figure(13,facecolor='white')
plt.savefig('3D-12hr.png')
plt.figure(14,facecolor='white')
plt.savefig('3D-24hr.png')
plt.figure(15,facecolor='white')
plt.savefig('3D-36hr.png')
plt.figure(16,facecolor='white')
plt.savefig('3D-48hr.png')
plt.figure(17,facecolor='white')
plt.savefig('3D-60hr.png')
plt.figure(18,facecolor='white')
plt.savefig('3D-72hr.png')
plt.figure(19,facecolor='white')
plt.savefig('3D-84hr.png')
plt.figure(20,facecolor='white')
plt.savefig('3D-96hr.png')
plt.figure(21,facecolor='white')
plt.savefig('3D-108hr.png')
plt.figure(22,facecolor='white')
plt.savefig('3D-120hr.png')
#plt.show()
Modules
The following are modules that were imported into the main program, main.py.
advectionx.py
This module decided the U wind velocity at each particle.
#advection-x
import numpy as np
from parameters import *
import scipy.interpolate as interp
R_advection = np.zeros((num_particles))
def advectionx(R,A,T,t):
for i in xrange(0,num_particles):
if A[i] == 1 and R[i,0] > 0.0:
alt_part = R[i,0]
lat_part = R[i,1]
lon_part = R[i,2]
a = int(R[i,0]/29000.*29.)
b = int((R[i,1]+90.)/180.*73.)
c = int(R[i,2]/360.*144.)
if a == 29:
a_high = a
else:
a_high = a+1
alt_low = alt_coor[a]
alt_high = alt_coor[a_high]
lat_down = lat_coor[b]
lat_up = lat_coor[b+1]
lon_left = lon_coor[c]
lon_right = lon_coor[c+1]
#trilinear interpolation
#time 1
interp_x1 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T,a,b,c+1] - U[T,a,b,c]) + U[T,a,b,c]
interp_x2 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T,a,b+1,c+1] - U[T,a,b+1,c]) + U[T,a,b+1,c]
interp_x3 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T,a_high,b,c+1] - U[T,a_high,b,c]) + U[T,a_high,b,c]
interp_x4 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T,a_high,b+1,c+1] - U[T,a_high,b+1,c]) + U[T,a_high,b+1,c]
interp_y1 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x2 - interp_x1) + interp_x1
interp_y2 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x4 - interp_x3) + interp_x3
if a_high == a:
U_interp_1 = (interp_y1+interp_y2)/2.0
else:
U_interp_1 = (alt_part - alt_low) / (alt_high - alt_low) * (interp_y2 - interp_y1) + interp_y1
#time 2
interp_x1 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T+1,a,b,c+1] - U[T+1,a,b,c]) + U[T+1,a,b,c]
interp_x2 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T+1,a,b+1,c+1] - U[T+1,a,b+1,c]) + U[T+1,a,b+1,c]
interp_x3 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T+1,a_high,b,c+1] - U[T+1,a_high,b,c]) + U[T+1,a_high,b,c]
interp_x4 = (lon_part - lon_left) / (lon_right - lon_left) * (U[T+1,a_high,b+1,c+1] - U[T+1,a_high,b+1,c]) + U[T+1,a_high,b+1,c]
interp_y1 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x2 - interp_x1) + interp_x1
interp_y2 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x4 - interp_x3) + interp_x3
if a_high == a:
U_interp_2 = (interp_y1+interp_y2)/2.0
else:
U_interp_2 = (alt_part - alt_low) / (alt_high - alt_low) * (interp_y2 - interp_y1) + interp_y1
R_advection[i] = (U_interp_2 - U_interp_1) / 1.0 * (t - float(T)) + U_interp_1
else:
R_advection[i] = 0.0
return R_advection
advectiony.py
This module decided the V wind velocity at each particle.
#advection
import numpy as np
from parameters import *
R_advection = np.zeros((num_particles))
def advectiony(R,A,T,t):
for i in xrange(0,num_particles):
if A[i] == 1 and R[i,0] > 0.0:
alt_part = R[i,0]
lat_part = R[i,1]
lon_part = R[i,2]
a = int(R[i,0]/29000.*29.)
b = int((R[i,1]+90.)/180.*73.)
c = int(R[i,2]/360.*144.)
if a == 29:
a_high = a
else:
a_high = a+1
alt_low = alt_coor[a]
alt_high = alt_coor[a_high]
lat_down = lat_coor[b]
lat_up = lat_coor[b+1]
lon_left = lon_coor
lon_right = lon_coor
#time 1
interp_x1 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T,a,b,c+1] – V[T,a,b,c]) + V[T,a,b,c]
interp_x2 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T,a,b+1,c+1] – V[T,a,b+1,c]) + V[T,a,b+1,c]
interp_x3 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T,a_high,b,c+1] – V[T,a_high,b,c]) + V[T,a_high,b,c]
interp_x4 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T,a_high,b+1,c+1] – V[T,a_high,b+1,c]) + V[T,a_high,b+1,c]
interp_y1 = (lat_part – lat_down) / (lat_up – lat_down) * (interp_x2 – interp_x1) + interp_x1
interp_y2 = (lat_part – lat_down) / (lat_up – lat_down) * (interp_x4 – interp_x3) + interp_x3
if a_high == a:
V_interp_1 = (interp_y1+interp_y2)/2.0
else:
V_interp_1 = (alt_part – alt_low) / (alt_high – alt_low) * (interp_y2 – interp_y1) + interp_y1
#time 2
interp_x1 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T+1,a,b,c+1] – V[T+1,a,b,c]) + V[T+1,a,b,c]
interp_x2 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T+1,a,b+1,c+1] – V[T+1,a,b+1,c]) + V[T+1,a,b+1,c]
interp_x3 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T+1,a_high,b,c+1] – V[T+1,a_high,b,c]) + V[T+1,a_high,b,c]
interp_x4 = (lon_part – lon_left) / (lon_right – lon_left) * (V[T+1,a_high,b+1,c+1] – V[T+1,a_high,b+1,c]) + V[T+1,a_high,b+1,c]
interp_y1 = (lat_part – lat_down) / (lat_up – lat_down) * (interp_x2 – interp_x1) + interp_x1
interp_y2 = (lat_part – lat_down) / (lat_up – lat_down) * (interp_x4 – interp_x3) + interp_x3
if a_high == a:
V_interp_2 = (interp_y1+interp_y2)/2.0
else:
V_interp_2 = (alt_part – alt_low) / (alt_high – alt_low) * (interp_y2 – interp_y1) + interp_y1
R_advection[i] = (V_interp_2 – V_interp_1) / 1.0 * (t – float(T)) + V_interp_1
else:
R_advection[i] = 0.0
return R_advection
advectionz.py
This module decided the W wind velocity at each particle.
#advection
import numpy as np
from parameters import *
R_advection = np.zeros((num_particles))
def advectionz(R,A,T,t):
for i in xrange(0,num_particles):
if A[i] == 1 and R[i,0] > 0.0:
alt_part = R[i,0]
lat_part = R[i,1]
lon_part = R[i,2]
a = int(R[i,0]/29000.*29.)
b = int((R[i,1]+90.)/180.*73.)
c = int(R[i,2]/360.*144.)
if a == 29:
a_high = a
else:
a_high = a+1
alt_low = alt_coor[a]
alt_high = alt_coor[a_high]
lat_down = lat_coor[b]
lat_up = lat_coor[b+1]
lon_left = lon_coor[c]
lon_right = lon_coor[c+1]
#trilinear interpolation
#time 1
interp_x1 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T,a,b,c+1] - W[T,a,b,c]) + W[T,a,b,c]
interp_x2 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T,a,b+1,c+1] - W[T,a,b+1,c]) + W[T,a,b+1,c]
interp_x3 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T,a_high,b,c+1] - W[T,a_high,b,c]) + W[T,a_high,b,c]
interp_x4 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T,a_high,b+1,c+1] - W[T,a_high,b+1,c]) + W[T,a_high,b+1,c]
interp_y1 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x2 - interp_x1) + interp_x1
interp_y2 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x4 - interp_x3) + interp_x3
if a_high == a:
W_interp_1 = (interp_y1+interp_y2)/2.0
else:
W_interp_1 = (alt_part - alt_low) / (alt_high - alt_low) * (interp_y2 - interp_y1) + interp_y1
#time 2
interp_x1 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T+1,a,b,c+1] - W[T+1,a,b,c]) + W[T+1,a,b,c]
interp_x2 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T+1,a,b+1,c+1] - W[T+1,a,b+1,c]) + W[T+1,a,b+1,c]
interp_x3 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T+1,a_high,b,c+1] - W[T+1,a_high,b,c]) + W[T+1,a_high,b,c]
interp_x4 = (lon_part - lon_left) / (lon_right - lon_left) * (W[T+1,a_high,b+1,c+1] - W[T+1,a_high,b+1,c]) + W[T+1,a_high,b+1,c]
interp_y1 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x2 - interp_x1) + interp_x1
interp_y2 = (lat_part - lat_down) / (lat_up - lat_down) * (interp_x4 - interp_x3) + interp_x3
if a_high == a:
W_interp_2 = (interp_y1+interp_y2)/2.0
else:
W_interp_2 = (alt_part - alt_low) / (alt_high - alt_low) * (interp_y2 - interp_y1) + interp_y1
R_advection[i] = (W_interp_2 - W_interp_1) / 1.0 * (t - float(T)) + W_interp_1
else:
R_advection[i] = 0.0
return R_advection
degree_conversion.py
This module converted a velocity in meters per second to a velocity in decimal degrees per second.
#settling
import numpy as np
from parameters import *
a = 6378137.0
def degree_to_dx(R):
h = R[:,0]
lat = R[:,1]
lon = R[:,2]
dx_lat = np.pi / 180. * (a + h)
dx_lon = np.pi / 180. * (a + h) * np.cos(lat * np.pi / 180. )
return dx_lat, dx_lon
dispersionx.py
This module simulates the dispersion in the x-direction.
#dispersion
import numpy as np
from parameters import *
R_dispersion = np.zeros((num_particles))
def dispersionx(A):
c_h = np.sqrt(2.0 * K_h /dt)
for i in xrange(0,num_particles):
if A[i] == 1:
R_dispersion[i] = np.random.normal(0.0,c_h)
else:
R_dispersion[i] = 0.0
return R_dispersion
dispersiony.py
This module simulates the dispersion in the y-direction.
#dispersion
import numpy as np
from parameters import *
R_dispersion = np.zeros((num_particles))
def dispersiony(A):
c_h = np.sqrt(2.0 * K_h /dt)
for i in xrange(0,num_particles):
if A[i] == 1:
R_dispersion[i] = np.random.normal(0.0,c_h)
else:
R_dispersion[i] = 0.0
return R_dispersion
dispersionz.py
This module simulates the dispersion in the z-direction.
#dispersion
import numpy as np
from parameters import *
R_dispersion = np.zeros((num_particles))
def dispersionz(A):
c_v = np.sqrt(2.0 * K_v /dt)
for i in xrange(0,num_particles):
if A[i] == 1:
R_dispersion[i] = np.random.normal(0.0,c_v)
else:
R_dispersion[i] = 0.0
return R_dispersion
initial.py
This module sets up the initial position of the particles.
import numpy as np
from parameters import *
R = np.zeros((num_particles,3))
def initial(sigma_lat,sigma_lon):
for i in xrange(0,num_particles):
R[i,0] = 15000.* np.random.uniform()
R[i,1] = np.random.normal(start_lat,sigma_lat)
R[i,2] = np.random.normal(start_lon,sigma_lon)
return R
landed.py
This module determines whether the particle is in the atmosphere or on the ground.
#landed
import numpy as np
from parameters import *
def landed(H,A):
for i in xrange(0,num_particles):
if A[i] == 1.0 and H[i] <= 0.0:
A[i] = 0.0
H[i] = 0.0
if H[i] > alt_coor[29]:
H[i] = alt_coor[29]
if H[i] <= 0.0:
H[i] = 0.0
return A, H
parameters.py
This module sets up our parameters.
import numpy as np
#parameters
num_particles = 1000
time = 5.0 #days
dt = 10.0 * 60.#minutes
time_steps = int(time * 24.0 * 3600. / dt) + 1
D_log = np.zeros((num_particles))
for i in xrange(0,num_particles):
D_log[i] = np.random.normal(-5.0,2.5)
D = 10. ** (D_log)
K_h = 2.0 * (10. **4.)
K_v = 10.
start_lat = 15.141667
start_lon = 120.35
sigma_lon = .01
sigma_lat = .01
U = np.load('u_wind.npy')
V = np.load('v_wind.npy')
W = np.load('w_wind.npy')
alt_coor = np.linspace(0,29000,30)
lat_coor = np.linspace(-90.,90.,73)
lon_coor = np.linspace(0.,360.,144)
settling.py
This module estimates the settling of the particles as a function of its diameter.
#settling
import numpy as np
from parameters import *
R_settling = np.zeros((num_particles))
constant = 1.09 * (10. ** 9.)
def settling(D,A):
for i in xrange(0,num_particles):
if A[i] == 1:
R_settling[i] = 2.0 / 9.0 * constant * (D[i] ** 2.0)
else:
R_settling[i] = 0.0
return R_settling
variables.py
This module sets up numpy arrays for our variables.
#variables import numpy as np from parameters import * #particle #, location of paritcles Rnew = np.zeros((num_particles,3)) #lat,lon,alt Rold = np.zeros((num_particles,3)) #lat,lon,alt R_step1 = np.zeros((num_particles,3)) #lat,lon,alt R_step2 = np.zeros((num_particles,3)) #lat,lon,alt R_step3 = np.zeros((num_particles,3)) #lat,lon,alt R_step4 = np.zeros((num_particles,3)) #lat,lon,alt Active = np.ones((num_particles)) #lat,lon,alt