Update toy_model.py

Corrected geometric mistake with Coriolis, tidied up optional oceans, changed the order of operation such that all thermal processes take place before updating the pressure

toy_model.py

@@ -9,7 +9,7 @@ import time, sys
 
 # define temporal parameters, including the length of time between calculation of fields and the length of a day on the planet (used for calculating Coriolis as well)
 day = 60*60*24
-dt = 60*1														###### <----- TIMESTEP
+dt = 60*6														###### <----- TIMESTEP
 
 # power incident on (lat,lon) at time t
 def solar(insolation, lat, lon, t):
@@ -41,13 +41,15 @@ u = np.zeros((nlat,nlon))
 v = np.zeros((nlat,nlon))
 air_density = np.zeros_like(air_pressure) + 1.3
 
-# if including an ocean, uncomment the below
+# custom oceans with lower albedo and higher heat capacity
-# albedo[5:55,9:20] = 0.2
+ocean = False
-# albedo[23:50,45:70] = 0.2
+if ocean:
-# albedo[2:30,85:110] = 0.2
+	albedo[5:55,9:20] = 0.2
-# heat_capacity_earth[5:55,9:20] = 1E6
+	albedo[23:50,45:70] = 0.2
-# heat_capacity_earth[23:50,45:70] = 1E6
+	albedo[2:30,85:110] = 0.2
-# heat_capacity_earth[2:30,85:110] = 1E6
+	heat_capacity_earth[5:55,9:20] = 1E6
+	heat_capacity_earth[23:50,45:70] = 1E6
+	heat_capacity_earth[2:30,85:110] = 1E6
 
 # define physical constants
 epsilon = 0.75
@@ -71,12 +73,12 @@ coriolis = np.zeros(nlat)	# also define the coriolis parameter here
 angular_speed = 2*np.pi/day
 for i in range(nlat):
 	dx[i] = dy*np.cos(lat[i]*np.pi/180)
-	coriolis[i] = day*np.cos(lat[i]*np.pi/180)
+	coriolis[i] = angular_speed*np.sin(lat[i]*np.pi/180)
 
 # define various useful differential functions:
 # gradient of scalar field a in the local x direction at point i,j
 def scalar_gradient_x(a,i,j):
-	return 0#(a[i,(j+1)%nlon]-a[i,(j-1)%nlon])/dx[i]
+	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
 def scalar_gradient_y(a,i,j):
 	if i == 0 or i == nlat-1:
@@ -112,14 +114,14 @@ ax[1].set_title('Atmosphere temperature')
 plt.ion()
 plt.show()
 
-# if you want to include advection set this to be True (currently this breaks the model!)
+# if you want to include advection set this to be True
 advection = True
 
 while True:
 
 	# print current time in simulation to command line
-	print("t = " + str(round(24*t/day,2)) + " days", end='\r')
+	print("t = " + str(round(t/day,2)) + " days", end='\r')
-	print(u.max())
+	# print(u.max(),air_density.max(),air_density.min())
 
 	# calculate change in temperature of ground and atmosphere due to radiative imbalance
 	for i in range(nlat):
@@ -127,6 +129,19 @@ while True:
 			temperature_planet[i,j] += dt*(albedo[i,j]*solar(insolation,lat[i],lon[j],t) + epsilon*sigma*temperature_atmosp[i,j]**4 - sigma*temperature_planet[i,j]**4)/heat_capacity_earth[i,j]
 			temperature_atmosp[i,j] += dt*(epsilon*sigma*temperature_planet[i,j]**4 - 2*epsilon*sigma*temperature_atmosp[i,j]**4)/heat_capacity_atmos
 
+	if advection:
+		boundary = 10
+		# allow for thermal advection in the atmosphere, and heat diffusion in the atmosphere and the ground
+		atmosp_addition = dt*divergence_with_scalar(temperature_atmosp)
+		temperature_atmosp[boundary:-boundary,:] -= atmosp_addition[boundary:-boundary,:]
+
+		# as density is now variable, allow for mass advection
+		density_addition = dt*divergence_with_scalar(air_density)
+		air_density[boundary:-boundary,:boundary] -= density_addition[boundary:-boundary,:boundary]
+
+	temperature_atmosp += dt*(thermal_diffusivity_air*laplacian(temperature_atmosp))
+	temperature_planet += dt*(thermal_diffusivity_roc*laplacian(temperature_planet))
+	
 	# update air pressure
 	air_pressure = air_density*specific_gas*temperature_atmosp
 
@@ -136,23 +151,12 @@ while True:
 	# calculate acceleration of atmosphere using primitive equations on beta-plane
 	for i in np.arange(1,nlat-1):
 		for j in range(nlon):
-			u_temp[i,j] += 0.001*dt*( -u[i,j]*scalar_gradient_x(u,i,j) - v[i,j]*scalar_gradient_y(u,i,j) + coriolis[i]*v[i,j] - scalar_gradient_x(air_pressure,i,j)/air_density[i,j] )
+			u_temp[i,j] += dt*(  - scalar_gradient_x(air_pressure,i,j)/air_density[i,j] + coriolis[i]*v[i,j]  - u[i,j]*scalar_gradient_x(u,i,j) - v[i,j]*scalar_gradient_y(u,i,j))
-			v_temp[i,j] += 0.001*dt*( -u[i,j]*scalar_gradient_x(v,i,j) - v[i,j]*scalar_gradient_y(v,i,j) - coriolis[i]*u[i,j] - scalar_gradient_y(air_pressure,i,j)/air_density[i,j] )
+			v_temp[i,j] += dt*(  - scalar_gradient_y(air_pressure,i,j)/air_density[i,j] - coriolis[i]*u[i,j]  - u[i,j]*scalar_gradient_x(v,i,j) - v[i,j]*scalar_gradient_y(v,i,j))
 
 	u += u_temp
 	v += v_temp
 
-	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_atmosp) + divergence_with_scalar(temperature_atmosp))
-		# temperature_atmosp[5:-5,:] += atmosp_addition[5:-5,:]
-
-		# as density is now variable, allow for mass advection
-		density_addition = dt*divergence_with_scalar(air_density)
-		# air_density[5:-5,:5] += density_addition[5:-5,:5]
-
-		temperature_planet += dt*(thermal_diffusivity_roc*laplacian(temperature_planet))
-	
 	# update plot
 	test = ax[0].contourf(lon_plot, lat_plot, temperature_planet, cmap='seismic')
 	ax[1].contourf(lon_plot, lat_plot, temperature_atmosp, cmap='seismic')