Smoothing added

Added FFT function to smooth velocity and temperature fields.
Also minor graphical improvements.
Added surface friction
Clamped velocities over pole

claude_low_level_library.pyx

@@ -13,7 +13,7 @@ cdef float inv_90 = np.pi/90
 def scalar_gradient_x(np.ndarray a,np.ndarray dx,np.int_t nlon,np.int_t i,np.int_t j,np.int_t k):
 	return (a[i,(j+1)%nlon,k]-a[i,(j-1)%nlon,k])/dx[i]
 
-def scalar_gradient_x_2D(a,dx,nlon,i,j):
+def scalar_gradient_x_2D(np.ndarray a,np.ndarray dx,np.int_t nlon,np.int_t i,np.int_t j):
 	return (a[i,(j+1)%nlon]-a[i,(j-1)%nlon])/dx[i]
 
 # gradient of scalar field a in the local y direction at point i,j
@@ -25,7 +25,7 @@ def scalar_gradient_y(np.ndarray a,DTYPE_f dy,np.int_t nlat,np.int_t i,np.int_t
 	else:
 		return (a[i+1,j,k]-a[i-1,j,k])/dy
 
-def scalar_gradient_y_2d(a,dy,nlat,i,j):
+def scalar_gradient_y_2D(np.ndarray a,DTYPE_f dy,np.int_t nlat,np.int_t i,np.int_t j):
 	if i == 0:
 		return 2*(a[i+1,j]-a[i,j])/dy
 	elif i == nlat-1:
@@ -33,7 +33,7 @@ def scalar_gradient_y_2d(a,dy,nlat,i,j):
 	else:
 		return (a[i+1,j]-a[i-1,j])/dy
 
-def scalar_gradient_z(a,dz,i,j,k):
+def scalar_gradient_z(np.ndarray a,np.ndarray dz,np.int_t i,np.int_t j,np.int_t k):
 	cdef np.int_t nlevels = len(dz)
 	if k == 0:
 		return (a[i,j,k+1]-a[i,j,k])/dz[k]
@@ -52,7 +52,8 @@ def scalar_gradient_z_1D(np.ndarray a,np.ndarray dz,np.int_t k):
 		return (a[k+1]-a[k-1])/(2*dz[k])
 
 def surface_optical_depth(DTYPE_f lat):
-	return 4 + np.cos(lat*inv_90)*2.5
+	cdef DTYPE_f inv_90
+	return 4 + np.cos(lat*inv_90)*2
 
 def thermal_radiation(DTYPE_f a):
 	cdef DTYPE_f sigma = 5.67E-8
@@ -83,5 +84,5 @@ def solar(DTYPE_f insolation,DTYPE_f  lat,DTYPE_f lon,np.int_t t,DTYPE_f  day,DT
 		else:
 			return value
 
-def profile(a):
+def profile(np.ndarray a):
 	return np.mean(np.mean(a,axis=0),axis=0)

claude_top_level_library.pyx

@@ -8,32 +8,52 @@ cimport cython
 ctypedef np.float64_t DTYPE_f
 
 # laplacian of scalar field a
-# def laplacian(a):
+def laplacian_2D(np.ndarray a,np.ndarray dx,DTYPE_f dy):
-# 	output = np.zeros_like(a)
+	cdef np.ndarray output = np.zeros_like(a)
-# 	if output.ndim == 2:
+	cdef np.int_t nlat,nlon,i,j
-# 		for i in np.arange(1,nlat-1):
+	cdef DTYPE_f inv_dx, inv_dy
-# 			for j in range(nlon):
+	nlat = a.shape[0]
-# 				output[i,j] = (scalar_gradient_x_2D(a,i,(j+1)%nlon) - scalar_gradient_x_2D(a,i,(j-1)%nlon))/dx[i] + (scalar_gradient_y_2D(a,i+1,j) - scalar_gradient_y_2D(a,i-1,j))/dy
+	nlon = a.shape[1]
-# 		return output
+	inv_dy = 1/dy
-# 	if output.ndim == 3:
+	for i in np.arange(1,nlat-1):
-# 		for i in np.arange(1,nlat-1):
+		inv_dx = dx[i]
-# 			for j in range(nlon):
+		for j in range(nlon):
-# 				for k in range(nlevels-1):
+			output[i,j] = (low_level.scalar_gradient_x_2D(a,dx,nlon,i,j) - low_level.scalar_gradient_x_2D(a,dx,nlon,i,j))*inv_dx + (low_level.scalar_gradient_y_2D(a,dy,nlat,i+1,j) - low_level.scalar_gradient_y_2D(a,dy,nlat,i-1,j))*inv_dy
-# 					output[i,j,k] = (scalar_gradient_x(a,i,(j+1)%nlon,k) - scalar_gradient_x(a,i,(j-1)%nlon,k))/dx[i] + (scalar_gradient_y(a,i+1,j,k) - scalar_gradient_y(a,i-1,j,k))/dy + (scalar_gradient_z(a,i,j,k+1)-scalar_gradient_z(a,i,j,k-1))/(2*dz[k])
+	return output
-# 		return output
+
+def laplacian_3D(np.ndarray a,np.ndarray dx,DTYPE_f dy,np.ndarray dz):
+	cdef np.ndarray output = np.zeros_like(a)
+	cdef np.int_t nlat,nlon,nlevels,i,j,k
+	cdef DTYPE_f inv_dx, inv_dy
+	nlat = a.shape[0]
+	nlon = a.shape[1]
+	nlevels = a.shape[2]
+	inv_dy = 1/dy
+	for i in np.arange(1,nlat-1):
+		inv_dx = 1/dx[i]
+		for j in range(nlon):
+			for k in range(nlevels-1):
+				output[i,j,k] = (low_level.scalar_gradient_x(a,dx,nlon,i,j,k) - low_level.scalar_gradient_x(a,dx,nlon,i,j,k))*inv_dx + (low_level.scalar_gradient_y(a,dy,nlat,i+1,j,k) - low_level.scalar_gradient_y(a,dy,nlat,i-1,j,k))*inv_dy + (low_level.scalar_gradient_z(a,dz,i,j,k+1)-low_level.scalar_gradient_z(a,dz,i,j,k-1))/(2*dz[k])
+	return output
 
 # divergence of (a*u) where a is a scalar field and u is the atmospheric velocity field
-def divergence_with_scalar(a,u,v,dx,dy):
+def divergence_with_scalar(np.ndarray a,np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray dx,DTYPE_f dy,np.ndarray dz):
-	output = np.zeros_like(a)
+	cdef np.ndarray output = np.zeros_like(a)
-	nlat, nlon, nlevels = output.shape[:]
+	cdef np.ndarray au, av, aw
+	cdef np.int_t nlat, nlon, nlevels, i, j, k 
+
+	nlat = output.shape[0]
+	nlon = output.shape[1]
+	nlevels = output.shape[2]
 
 	au = a*u
 	av = a*v
+	aw = a*w
 
 	for i in range(nlat):
 		for j in range(nlon):
 			for k in range(nlevels):
-				output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) #+ 0.1*scalar_gradient_z(a*w,i,j,k)
+				output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + 0.05*low_level.scalar_gradient_z(aw,dz,i,j,k)
 	return output
 
 def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_atmos, np.ndarray air_pressure, np.ndarray air_density, np.ndarray heat_capacity_earth, np.ndarray albedo, DTYPE_f insolation, np.ndarray lat, np.ndarray lon, np.ndarray heights, np.ndarray dz, np.int_t t, np.int_t dt, DTYPE_f day, DTYPE_f year, DTYPE_f axial_tilt):
@@ -41,6 +61,7 @@ def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_a
 	cdef np.int_t nlat,nlon,nlevels,i,j
 	cdef DTYPE_f fl = 0.1
 	cdef np.ndarray upward_radiation,downward_radiation,optical_depth,Q, pressure_profile, density_profile
+	cdef DTYPE_f sun_lat
 
 	nlat = temperature_atmos.shape[0]	
 	nlon = temperature_atmos.shape[1]
@@ -52,11 +73,12 @@ def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_a
 	Q = np.zeros(nlevels)
 
 	for i in range(nlat):
+		sun_lat = low_level.surface_optical_depth(lat[i])
 		for j in range(nlon):
 			# calculate optical depth
 			pressure_profile = air_pressure[i,j,:]
 			density_profile = air_density[i,j,:]
-			optical_depth = low_level.surface_optical_depth(lat[i])*(fl*(pressure_profile/pressure_profile[0]) + (1-fl)*(pressure_profile/pressure_profile[0])**4)
+			optical_depth = sun_lat*(fl*(pressure_profile/pressure_profile[0]) + (1-fl)*(pressure_profile/pressure_profile[0])**4)
 			
 			# calculate upward longwave flux, bc is thermal radiation at surface
 			upward_radiation[0] = low_level.thermal_radiation(temperature_world[i,j])
@@ -83,21 +105,28 @@ def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_a
 	
 	return temperature_world, temperature_atmos
 
-def velocity_calculation(u,v,air_pressure,old_pressure,air_density,coriolis,gravity,dx,dy,dt):
+def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray air_pressure,np.ndarray old_pressure,np.ndarray air_density,np.ndarray coriolis,DTYPE_f gravity,np.ndarray dx,DTYPE_f dy,DTYPE_f dt):
 	# introduce temporary arrays to update velocity in the atmosphere
-	u_temp = np.zeros_like(u)
+	cdef np.ndarray u_temp = np.zeros_like(u)
-	v_temp = np.zeros_like(v)
+	cdef np.ndarray v_temp = np.zeros_like(v)
-	w_temp = np.zeros_like(u)
+	cdef np.ndarray w_temp = np.zeros_like(u)
 
-	nlat,nlon,nlevels = air_pressure.shape[:]
+	cdef np.int_t nlat,nlon,nlevels,i,j,k
+	cdef DTYPE_f inv_density,inv_dt_grav
+
+	nlat = air_pressure.shape[0]
+	nlon = air_pressure.shape[1]
+	nlevels = air_pressure.shape[2]
+	inv_dt_grav = 1/(dt*gravity)
 
 	# calculate acceleration of atmosphere using primitive equations on beta-plane
-	for i in np.arange(1,nlat-1):
+	for i in np.arange(3,nlat-3).tolist():
 		for j in range(nlon):
 			for k in range(nlevels):
-				u_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(u,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(u,dy,nlat,i,j,k) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)/air_density[i,j,k] )
+				inv_density = 1/air_density[i,j,k]
-				v_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(v,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(v,dy,nlat,i,j,k) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(air_pressure,dy,nlat,i,j,k)/air_density[i,j,k] )
+				u_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(u,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(u,dy,nlat,i,j,k) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)*inv_density )
-				w_temp[i,j,k] += -(air_pressure[i,j,k]-old_pressure[i,j,k])/(dt*air_density[i,j,k]*gravity)
+				v_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(v,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(v,dy,nlat,i,j,k) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(air_pressure,dy,nlat,i,j,k)*inv_density )
+				w_temp[i,j,k] += -(air_pressure[i,j,k]-old_pressure[i,j,k])*inv_density*inv_dt_grav
 
 	u += u_temp
 	v += v_temp
@@ -107,4 +136,18 @@ def velocity_calculation(u,v,air_pressure,old_pressure,air_density,coriolis,grav
 	u *= 0.95
 	v *= 0.95
 
-	return u,v,w
+	u[:,:,0] *= 0.8
+	v[:,:,0] *= 0.8
+
+	return u,v,w
+
+def smoothing_3D(np.ndarray a,DTYPE_f smooth_parameter):
+	nlat = a.shape[0]
+	nlon = a.shape[1]
+	nlevels = a.shape[2]
+	smooth_parameter *= 0.5
+	test = np.fft.fftn(a)
+	test[int(nlat*smooth_parameter):int(nlat*(1-smooth_parameter)),:,:] = 0
+	test[:,int(nlon*smooth_parameter):int(nlon*(1-smooth_parameter)),:] = 0
+	# test[:,:int(nlevels*smooth_parameter):int(nlevels*(1-smooth_parameter))] = 0
+	return np.fft.ifftn(test).real

toy_model.py

@@ -12,30 +12,34 @@ import claude_top_level_library as top_level
 ######## CONTROL ########
 
 day = 60*60*24					# define length of day (used for calculating Coriolis as well) (s)
-resolution = 3					# how many degrees between latitude and longitude gridpoints
+resolution = 5					# how many degrees between latitude and longitude gridpoints
-nlevels = 10					# how many vertical layers in the atmosphere
+nlevels = 20					# how many vertical layers in the atmosphere
-top = 50E3						# top of atmosphere (m)
+top = 40E3						# top of atmosphere (m)
 planet_radius = 6.4E6			# define the planet's radius (m)
 insolation = 1370				# TOA radiation from star (W m^-2)
 gravity = 9.81 					# define surface gravity for planet (m s^-2)
-axial_tilt = -23.5				# tilt of rotational axis w.r.t. solar plane
+axial_tilt = -23.5/2				# tilt of rotational axis w.r.t. solar plane
 year = 365*day					# length of year (s)
 
 dt_spinup = 60*137
-dt_main = 60*9
+dt_main = 60*7
-spinup_length = 5*day
+spinup_length = 3*day
 
 ###
 
 advection = True 				# if you want to include advection set this to be True
-advection_boundary = 8			# how many gridpoints away from poles to apply advection
+advection_boundary = 3			# how many gridpoints away from poles to apply advection
+smoothing_parameter_t = 0.9
+smoothing_parameter_u = 0.8
+smoothing_parameter_v = 0.8
+smoothing_parameter_w = 0.3
 
 save = True 					# save current state to file?
 load = True  					# load initial state from file?
 
 plot = False 					# display plots of output?
-level_plots = True 				# display plots of output on vertical levels?
+level_plots = False 			# display plots of output on vertical levels?
-nplots = 5						# how many levels you want to see plots of (evenly distributed through column)
+nplots = 3						# how many levels you want to see plots of (evenly distributed through column)
 
 ###########################
 
@@ -98,7 +102,6 @@ for i in range(nlat):
 
 specific_gas = 287
 thermal_diffusivity_roc = 1.5E-6
-sigma = 5.67E-8
 
 air_pressure = air_density*specific_gas*temperature_atmos
 
@@ -121,17 +124,34 @@ dz[-1] = dz[-2]
 
 #################### SHOW TIME ####################
 
+# INITIATE TIME
+t = 0
+
+if load:
+	# load in previous save file
+	temperature_atmos,temperature_world,u,v,w,t,air_pressure,air_density,albedo = pickle.load(open("save_file.p","rb"))
+
 if plot:
 	# set up plot
 	f, ax = plt.subplots(2,figsize=(9,9))
 	f.canvas.set_window_title('CLAuDE')
-	ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic')
+	test = ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic')
-	ax[1].contourf(lon_plot, lat_plot, temperature_atmos[:,:,0], cmap='seismic')
+	ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=1)
+	ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic',levels=15)
+	ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8)
+	ax[1].streamplot(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1)),np.transpose(np.mean(10*w,axis=1)),color='black',density=0.75)
 	plt.subplots_adjust(left=0.1, right=0.75)
 	ax[0].set_title('Ground temperature')
+	ax[0].set_xlim(lon.min(),lon.max())
 	ax[1].set_title('Atmosphere temperature')
+	ax[1].set_xlim(lat.min(),lat.max())
+	ax[1].set_ylim(heights.min(),heights.max())
 	cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])
 
+	f.colorbar(test, cax=cbar_ax)
+	cbar_ax.set_title('Temperature (K)')
+	f.suptitle( 'Time ' + str(round(t/day,2)) + ' days' )
+
 	# allow for live updating as calculations take place
 
 	if level_plots:
@@ -150,13 +170,7 @@ if plot:
 
 	plt.ion()
 	plt.show()
-
+	plt.pause(2)
-# INITIATE TIME
-t = 0
-
-if load:
-	# load in previous save file
-	temperature_atmos,temperature_world,u,v,w,t,air_density,albedo = pickle.load(open("save_file.p","rb"))
 
 while True:
 
@@ -170,15 +184,14 @@ while True:
 		velocity = True
 
 	# print current time in simulation to command line
-	print("+++ t = " + str(round(t/day,2)) + " days +++", end='\r')
+	print("+++ t = " + str(round(t/day,2)) + " days +++")
 	print('T: ',round(temperature_world.max()-273.15,1),' - ',round(temperature_world.min()-273.15,1),' C')
 	print('U: ',round(u.max(),2),' - ',round(u.min(),2),' V: ',round(v.max(),2),' - ',round(v.min(),2),' W: ',round(w.max(),2),' - ',round(w.min(),2))
-	# print(profile(air_density))
-	# print(profile(air_pressure)/100)
 
-	before_radiation = time.time()
+	# before_radiation = time.time()
 	temperature_world, temperature_atmos = top_level.radiation_calculation(temperature_world, temperature_atmos, air_pressure, air_density, heat_capacity_earth, albedo, insolation, lat, lon, heights, dz, t, dt, day, year, axial_tilt)
-	time_taken = float(round(time.time() - before_radiation,3))
+	temperature_atmos = top_level.smoothing_3D(temperature_atmos,smoothing_parameter_t)
+	# time_taken = float(round(time.time() - before_radiation,3))
 	# print('Radiation: ',str(time_taken),'s')
 
 	# update air pressure
@@ -189,31 +202,49 @@ while True:
 
 		before_velocity = time.time()
 		u,v,w = top_level.velocity_calculation(u,v,air_pressure,old_pressure,air_density,coriolis,gravity,dx,dy,dt)
-		time_taken = float(round(time.time() - before_velocity,3))
+		u = top_level.smoothing_3D(u,smoothing_parameter_u)
+		# v = top_level.smoothing_3D(v,smoothing_parameter_v)
+		# w = top_level.smoothing_3D(w,smoothing_parameter_w)
+		
+		u[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
+		v[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
+		w[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
+
+		u[:advection_boundary,:,:] = 0
+		v[:advection_boundary,:,:] = 0
+		w[:advection_boundary,:,:] = 0
+
+		u[-advection_boundary:,:,:] = 0
+		v[-advection_boundary:,:,:] = 0
+		w[-advection_boundary:,:,:] = 0
+
+		# time_taken = float(round(time.time() - before_velocity,3))
 		# print('Velocity: ',str(time_taken),'s')
 
-		before_advection = time.time()
+		# before_advection = time.time()
 		if advection:
 			# allow for thermal advection in the atmosphere, and heat diffusion in the atmosphere and the ground
 			# atmosp_addition = dt*(thermal_diffusivity_air*laplacian(temperature_atmos))
 
-			atmosp_addition = dt*top_level.divergence_with_scalar(temperature_atmos,u,v,dx,dy)
+			atmosp_addition = dt*top_level.divergence_with_scalar(temperature_atmos,u,v,w,dx,dy,dz)
 			temperature_atmos[advection_boundary:-advection_boundary,:,:] -= atmosp_addition[advection_boundary:-advection_boundary,:,:]
 			temperature_atmos[advection_boundary-1,:,:] -= 0.5*atmosp_addition[advection_boundary-1,:,:]
 			temperature_atmos[-advection_boundary,:,:] -= 0.5*atmosp_addition[-advection_boundary,:,:]
 
 			# as density is now variable, allow for mass advection
-			# density_addition = dt*divergence_with_scalar(air_density)
+			density_addition = dt*top_level.divergence_with_scalar(air_density,u,v,w,dx,dy,dz)
-			# air_density[advection_boundary:-advection_boundary,:,:] -= density_addition[advection_boundary:-advection_boundary,:,:]
+			air_density[advection_boundary:-advection_boundary,:,:] -= density_addition[advection_boundary:-advection_boundary,:,:]
-			# air_density[(advection_boundary-1),:,:] -= 0.5*density_addition[advection_boundary-1,:,:]
+			air_density[(advection_boundary-1),:,:] -= 0.5*density_addition[advection_boundary-1,:,:]
-			# air_density[-advection_boundary,:,:] -= 0.5*density_addition[-advection_boundary,:,:]
+			air_density[-advection_boundary,:,:] -= 0.5*density_addition[-advection_boundary,:,:]
+
+			# temperature_world -= dt*(thermal_diffusivity_roc*top_level.laplacian_2D(temperature_world,dx,dy))
 
-			# temperature_world += dt*(thermal_diffusivity_roc*laplacian(temperature_world))
 			
-			time_taken = float(round(time.time() - before_advection,3))
+			# time_taken = float(round(time.time() - before_advection,3))
 			# print('Advection: ',str(time_taken),'s')
+
 	
-	before_plot = time.time()
+	# before_plot = time.time()
 	if plot:
 
 		tropopause_height = np.zeros(nlat)
@@ -226,10 +257,11 @@ while True:
 			else:
 				while low_level.scalar_gradient_z_1D(np.mean(temperature_atmos[i,:,:],axis=0),dz,k) < 0: 
 					k += 1
+			tropopause_height[i] = heights[k]
 
 		# update plot
-		test = ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic')
+		test = ax[0].contourf(lon_plot, lat_plot, temperature_world, cmap='seismic',levels=15)
-		ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=1)
+		if velocity:	ax[0].streamplot(lon_plot, lat_plot, u[:,:,0], v[:,:,0], color='white',density=1)
 		ax[0].set_title('$\it{Ground} \quad \it{temperature}$')
 		ax[0].set_xlim((lon.min(),lon.max()))
 		ax[0].set_ylim((lat.min(),lat.max()))
@@ -237,10 +269,11 @@ while True:
 		ax[0].axhline(y=0,color='black',alpha=0.3)
 		ax[0].set_xlabel('Longitude')
 
-		ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic')
+		ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic',levels=15)
 		ax[1].plot(lat_plot,tropopause_height,color='black',linestyle='--',linewidth=3,alpha=0.5)
-		ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8)
+		if velocity:
-		ax[1].streamplot(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1)),np.transpose(np.mean(10*w,axis=1)),color='black',density=0.75)
+			ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8)
+			ax[1].streamplot(heights_plot, lat_z_plot, np.transpose(np.mean(v,axis=1)),np.transpose(np.mean(10*w,axis=1)),color='black',density=0.75)
 		ax[1].set_title('$\it{Atmospheric} \quad \it{temperature}$')
 		ax[1].set_xlim((-90,90))
 		ax[1].set_ylim((0,heights.max()))
@@ -252,10 +285,12 @@ while True:
 		f.suptitle( 'Time ' + str(round(24*t/day,2)) + ' hours' )
 		
 		if level_plots:
+			quiver_padding = int(50/resolution)
+			skip=(slice(None,None,2),slice(None,None,2))
 			for k, z in zip(range(nplots), level_plots_levels):	
 				z += 1
-				bx[k].contourf(lon_plot, lat_plot, temperature_atmos[:,:,z], cmap='seismic')
+				bx[k].contourf(lon_plot, lat_plot, temperature_atmos[:,:,z], cmap='seismic',levels=15)
-				bx[k].streamplot(lon_plot, lat_plot, u[:,:,z], v[:,:,z], color='white',density=1)
+				bx[k].streamplot(lon_plot, lat_plot, u[:,:,z], v[:,:,z], color='white',density=1.5)
 				bx[k].set_title(str(round(heights[z]/1E3))+' km')
 				bx[k].set_ylabel('Latitude')
 				bx[k].set_xlim((lon.min(),lon.max()))
@@ -270,7 +305,7 @@ while True:
 		if level_plots:
 			for k in range(nplots):
 				bx[k].cla()
-	time_taken = float(round(time.time() - before_plot,3))
+	# time_taken = float(round(time.time() - before_plot,3))
 	# print('Plotting: ',str(time_taken),'s')
 
 	# advance time by one timestep
@@ -281,7 +316,7 @@ while True:
 	print('Time: ',str(time_taken),'s')
 
 	if save:
-		pickle.dump((temperature_atmos,temperature_world,u,v,w,t,air_density,albedo), open("save_file.p","wb"))
+		pickle.dump((temperature_atmos,temperature_world,u,v,w,t,air_pressure,air_density,albedo), open("save_file.p","wb"))
 
 	if np.isnan(u.max()):
 		sys.exit()