import numpy as np
import matplotlib.pyplot as plt
import time
#climate model
def climate_rate(t,n,parameters):
sec_to_yr = 3600. * 24. * 365.25
MLD = 50.0
Cv_water = 4184.
rho_water = 1000.
beta_surf = sec_to_yr / rho_water / Cv_water / MLD
alpha = 1.48
Cv_air = 700.
beta_atms = alpha * sec_to_yr / 10000. / Cv_air
cloud_albedo = .225 / .60
clouds_boolean = parameters[0]
greenhouse_boolean = parameters [6]
if clouds_boolean ==0:
R_atms = .225 # -
else:
R_atms = cloud_albedo * np.interp(n[0],cloud_data[:,0],cloud_data[:,1])
Abs_atms = parameters[1] # -
E_sol = parameters[2] # W/m2
CO2 = parameters[3] # ppm
KH2O = parameters[4] # -
#0 is surface, 1 is atmosphere
T_surf = n[0]
T_atms = n[1]
#extra parameter
R_surf = 0.08465 + 0.38 * np.exp(-0.006 * (T_surf - 260.)**2.0)
H2O = KH2O * (0.6 + 0.5 * np.exp(5420. * (T_surf - 288.03)\
/(T_surf * 288.03)))
ems_atms = 0.76 + 0.030 * np.sqrt(CO2 / 320.0) + 0.1 * H2O
LSH_flux = 104.0 * ((T_surf - T_atms)/20.55)**0.5
#Surface
surf_in_1 = E_sol * (1.-R_surf-R_atms-Abs_atms)
surf_in_2 = 5.67 * (10.0 ** - 8.0)\
*ems_atms * 1.25 * (T_atms ** 4.0)
surf_out = 5.67 * (10.0 ** - 8.0)\
* (T_surf ** 4.0)\
+ LSH_flux
#Atmosphere
atms_in_1 = Abs_atms * E_sol
if greenhouse_boolean == 0:
atms_in_2 = 5.67 * (10.0 ** - 8.0)\
+ LSH_flux
else:
atms_in_2 = 5.67 * (10.0 ** - 8.0)\
* ems_atms * (T_surf ** 4.0)\
+ LSH_flux
atms_out = 5.67 * (10.0 ** - 8.0)\
* (2.0 * ems_atms * T_atms ** 4.0)
rate1 = (surf_in_1 + surf_in_2 - surf_out) * beta_surf
rate2 = (atms_in_1 + atms_in_2 - atms_out) * beta_atms
return rate1,rate2
#4th order runga kutta
def runga_kutta_4(function,parameters,Ni,ti,tf,dt):
N_size = np.shape(Ni)
last_N = N_size[0]+1
steps = int((tf-ti)/dt)+1
N = np.zeros((steps,last_N))
N[0,0] = ti
N[0,1:last_N] = Ni
for i in xrange(1,steps):
k1 = dt * np.array(function(N[i-1,0],N[i-1,1:last_N],parameters))
k2 = dt * np.array(function(N[i-1,0]+dt/2.0,(N[i-1,1:last_N]+k1/2.0),parameters))
k3 = dt *np.array(function(N[i-1,0]+dt/2.0,(N[i-1,1:last_N]+k2/2.0),parameters))
k4 = dt * np.array(function(N[i-1,0]+dt,(N[i-1,1:last_N]+k3),parameters))
N[i,0] = N[i-1,0]+ dt
N[i,1:last_N] = N[i-1,1:last_N] + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0
return N
##parameters
##parameters[0] - Clouds on? # - 0 off, 1 on
##Abs_atms = parameters[1] # -
##E_sol = parameters[2] # W/m2
##CO2 = parameters[3] # ppm
##KH2O = parameters[4] # -
##parameters[5] #lag time
##parameters [6] GreenHouse on? # 0 off, 1 on
###RUN?
run_steady = 0
run_sensitivity = 0
run_response = 0
run_cloud = 0
run_greenhouse1 = 0
run_greenhouse2 = 0
#Steady State
if run_steady == 1:
soln1 = runga_kutta_4(climate_rate,[0 , 0.196, 342., 320.0, 1.0, -1.0, 1],[288.0667,267.4875],0.,30.,0.001)
plt.figure(1,facecolor='white')
plt.plot(soln1[:,0],soln1[:,1],label = 'surface')
plt.plot(soln1[:,0],soln1[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(0,30)
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('steady_state.png')
print 'steady state done'
#Sensitivity
if run_sensitivity == 1:
list_len = 51
final_T = np.zeros((list_len,2))
response = np.zeros((list_len,2))
E_sol_list = np.linspace(342.*.97,342.*1.03,list_len)
for i in xrange(0,list_len):
soln2 = runga_kutta_4(climate_rate,[0 , 0.196, E_sol_list[i], 320.0, 1.0, -1.0, 1],[288.0667,267.4875],0.,30.,0.001)
final_T[i,0] = soln2[-1,1]
final_T[i,1] = soln2[-1,2]
print str(int(float(i)/float(list_len)*100.)) + '%'
if i == 0:
plt.figure(2,facecolor='white')
plt.plot(soln2[:,0],soln2[:,1],label = 'surface')
plt.plot(soln2[:,0],soln2[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(0,30)
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('97%.png')
elif i == list_len-1:
plt.figure(3,facecolor='white')
plt.plot(soln2[:,0],soln2[:,1],label = 'surface')
plt.plot(soln2[:,0],soln2[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(0,30)
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('103%.png')
diff_surf = np.abs(soln2[0,1]-soln2[-1,1])
diff_atms = np.abs(soln2[0,2]-soln2[-1,2])
bingo1 = 0
bingo2 = 0
for t in xrange(0,30001):
if bingo1 == 0:
if np.abs(soln2[t,1]-soln2[0,1])/diff_surf > .95:
bingo1 = 1
response[i,0] = (float(t)*0.001)
if bingo2 == 0:
if np.abs(soln2[t,2]-soln2[0,2])/diff_atms > .95:
bingo2 = 1
response[i,1] = (float(t)*0.001)
plt.figure(4,facecolor='white')
plt.plot(E_sol_list,final_T[:,0],label = 'surface')
plt.plot(E_sol_list,final_T[:,1],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(E_sol_list[0],E_sol_list[-1])
plt.xlabel(r'Solar Input [$W/m^2$]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('sensitivity.png')
print 'sensitivity done'
plt.figure(5,facecolor='white')
plt.plot(E_sol_list,response[:,0],label = 'surface')
plt.plot(E_sol_list,response[:,1],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(E_sol_list[0],E_sol_list[-1])
plt.xlabel(r'Solar Input [$W/m^2$]')
plt.ylabel('response time [yr]')
plt.savefig('response.png')
if run_response == 1:
soln3 = runga_kutta_4(climate_rate,[0 , 0.196, 342., 320.0, 1.0, -1.0, 1],[288.0667,277.4875],0.,30.,0.001)
plt.figure(5,facecolor='white')
plt.plot(soln3[:,0],soln3[:,1],label = 'surface')
plt.plot(soln3[:,0],soln3[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.xlim(0,30)
plt.savefig('response2.png')
#Cloud
if run_cloud == 1:
cloud_data = np.zeros((7,2))
cloud_data[:,0] = (258., 268., 278., 288., 298., 308., 318.)
cloud_data[:,1] = (0.015,0.110,0.360,0.600,0.840,.920,.950)
soln4 = runga_kutta_4(climate_rate,[1 , 0.196, 342. * 1.03 , 320.0, 1.0, -1.0, 1],[288.0667,267.4875],0.,30.,0.001)
plt.figure(6,facecolor='white')
plt.plot(soln4[:,0],soln4[:,1],label = 'surface')
plt.plot(soln4[:,0],soln4[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(0,25)
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('cloud.png')
print 'clouds done'
#GreenHouse- turning it off
if run_greenhouse1 == 1:
soln5 = runga_kutta_4(climate_rate,[0 , 0.196, 342., 320.0, 1.0, -1.0, 0],[288.0667,267.4875],0.,30.,0.001)
plt.figure(7,facecolor='white')
plt.plot(soln5[:,0],soln5[:,1],label = 'surface')
plt.plot(soln5[:,0],soln5[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(0,30)
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('GreenHouse_off.png')
print 'GreenHouse1 done'
#GreenHouse- enhancing it
if run_greenhouse2 == 1:
soln6 = runga_kutta_4(climate_rate,[0 , 0.196, 342., 320.0 * 2.0 , 1.0, -1.0, 1],[288.0667,267.4875],0.,30.,0.001)
plt.figure(8,facecolor='white')
plt.plot(soln6[:,0],soln6[:,1],label = 'surface')
plt.plot(soln6[:,0],soln6[:,2],label = 'atmosphere')
plt.legend(loc='best')
plt.xlim(0,30)
plt.xlabel('Time [Years]')
plt.ylabel('Temperature [Celcius]')
plt.savefig('GreenHouse_enhance.png')
print 'GreenHouse2 done'
plt.show()