From df23337825e5aa5816b3db11af63d09ac76006ce Mon Sep 17 00:00:00 2001 From: Simon Clark Date: Mon, 15 Jun 2020 22:10:59 +0100 Subject: [PATCH] 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 | 56 ++++++++++++++++++++++++++++------------------------ 1 file changed, 30 insertions(+), 26 deletions(-) diff --git a/toy_model.py b/toy_model.py index d51d8c3..e29a8f3 100644 --- a/toy_model.py +++ b/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 -# albedo[5:55,9:20] = 0.2 -# albedo[23:50,45:70] = 0.2 -# albedo[2:30,85:110] = 0.2 -# heat_capacity_earth[5:55,9:20] = 1E6 -# heat_capacity_earth[23:50,45:70] = 1E6 -# heat_capacity_earth[2:30,85:110] = 1E6 +# custom oceans with lower albedo and higher heat capacity +ocean = False +if ocean: + albedo[5:55,9:20] = 0.2 + albedo[23:50,45:70] = 0.2 + albedo[2:30,85:110] = 0.2 + 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(u.max()) + print("t = " + str(round(t/day,2)) + " days", end='\r') + # 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] ) - 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] ) + 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] += 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')