File:Earth like nontilted ocean planet pco2 effect to ts 1 1 1 1.png

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3
2. 1. seasonal climlab energy balance model
3. python3/climblab code
5. 14.11.2023 0000.0004d2

import math import numpy as np import matplotlib.pyplot as plt from matplotlib import cm import climlab from climlab import constants as const from climlab.process.diagnostic import DiagnosticProcess from climlab.domain.field import Field, global_mean from scipy.interpolate import griddata from skimage.transform import resize

class tanalbedo(DiagnosticProcess):

   def __init__(self, **kwargs):
       super(tanalbedo, self).__init__(**kwargs)

       self.add_diagnostic('albedo')
       Ts = self.state['Ts']
       self._compute_fixed()

   def _compute_fixed(self):
       Ts = self.state['Ts']
       try:
           lon, lat = np.meshgrid(self.lon, self.lat)
       except:
           lat = self.lat
       phi = lat

       try:
           albedo=np.zeros(len(phi));
           albedo=0.42-0.20*np.tanh(0.052*(Ts-3))

       except:
           albedo = np.zeros_like(phi)

       dom = next(iter(self.domains.values()))
       self.albedo = Field(albedo, domain=dom)

   def _compute(self):
       self._compute_fixed()

       return {}

def run_ebm_00(Sk, co2): numyears=30 ##n no function here, run to equil numlat=18 numlev=6 plotvar=0 ## 1,2,3 lot temp, ice, mean albedo waterdepth=20 #S1=1365.2*1 au1=1.00 #Sk=1/math.pow(au1,2) ## relative sun constant to Earth now S1=1361.5*Sk #ecc=0.0167643, #long_peri=280.32687 #obliquity=23.459277 ecc=0 long_peri=0 obliquity=90 #co2=280 ##co2 amount ppmv #co2=280 diffu1=0.3 # meridional diffusivity in m**2/s albedo0=0.28 #orbit={'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 23.459277, 'S0':S1} orbit={'ecc': ecc, 'long_peri': long_peri, 'obliquity': obliquity, 'S0':S1} # creating EBM model #ebm= climlab.EBM(CO2=co2,orbit={'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 23.459277, 'S0':S1}) #ebm0= climlab.EBM_seasonal(water_depth=10.0, a0=0.3, num_lat=90, lum_lon=None, num_lev=10,num_lon=None #, orbit=orbit) ebm0= climlab.EBM_seasonal(water_depth=waterdepth, a0=albedo0, num_lat=numlat, lum_lon=None, num_lev=numlev,num_lon=None) ebm=climlab.process_like(ebm0) #ebm.step_forward() #print(ebm.diagnostics) #quit(-1) surface = ebm.domains['Ts'] # define new insolation and SW process ebm.remove_subprocess('insolation') insolation = climlab.radiation.DailyInsolation(domains=surface, orb = orbit, **ebm.param) insolation.S0=S1 ##sun = climlab.radiation.DailyInsolation(domains=model.Ts.domain) ebm.add_subprocess('insolation', insolation) #ebm.step_forward() #print(insolation.diagnostics) #print (insolation.insolation) #print (np.max(insolation.insolation)) ##print(insolation.S0) #quit(-1) ebm.remove_subprocess('albedo') alb = climlab.surface.albedo.StepFunctionAlbedo(state=ebm.state, Tf=-10, **ebm.param) #alb = climlab.surface.albedo.StepFunctionAlbedo(state=ebm.state, Tf=-20, **ebm.param) #alb = climlab.surface.albedo.ConstantAlbedo(domains=surface, **ebm.param) #alb = tanalbedo(state=ebm.state, **ebm.param) ebm.add_subprocess('albedo', alb) ebm.remove_subprocess('SW') SW = climlab.radiation.absorbed_shorwave.SimpleAbsorbedShortwave(insolation=insolation.insolation, state = ebm.state, albedo = alb.albedo, **ebm.param) ebm.add_subprocess('SW', SW) ebm.remove_subprocess('LW') LW = climlab.radiation.aplusbt.AplusBT_CO2(CO2=co2,state=ebm.state, **ebm.param) ebm.add_subprocess('LW', LW) #print(SW.diagnostics) #quit(-1) #ebm.CO2=co2 ebm.remove_subprocess('diffusion') D=diffu1 # meridional diffusivity in m**2/s #K = D / ebm.Tatm.domain.heat_capacity * const.a**2 K= D/ 700* const.a**2 diff = climlab.dynamics.MeridionalMoistDiffusion(state=ebm.state, timestep=ebm.timestep) ebm.add_subprocess('diffusion', diff) #print (ebm) ebm.step_forward() #ebm.diagnostics #ebm.integrate_years(numyears) #ebm.integrate_years(1) ebm.integrate_converge() #print(ebm.Ts) #plt.plot(ebm.Ts) #plt.show() num_steps_per_year = int(ebm.time['num_steps_per_year']) mean_year = np.empty(num_steps_per_year) for m in range(num_steps_per_year): ebm.step_forward() mean_year[m] = ebm.global_mean_temperature()

Tmean_year = np.mean(mean_year) print(round(Tmean_year,2)) return(Tmean_year)

1. co2s=[0,1,10,100,1000, 10000, 100000, 1e6]

1. co2s=[0,10,100, 10000, 1e6]

1. co2s=[0,1,3,10,30, 100,280, 1000, 3000, 10000, 30000, 100000, 300000, 1e6]

co2s=[0,0.01, 0.1, 0.3,1,3,10,30,50,100,280,400,600, 1000,1500, 2000,3000, 4000, 6000,10000, 30000, 100000, 300000, 1e6]

1. co2s=[0,1000, 1e6]

Tss=[] lenu=len(co2s)

for n in range(0,lenu): Ts=run_ebm_00(1.0, co2s[n]) print(co2s[n],Ts) Tss.append(Ts)

co2t=np.array(co2s) Tst=np.array(Tss)

plt.plot(np.log10(co2t), Tst, lw=3, color="#000000") plt.title("Earth-like ocean planet CO2 --- T", fontsize=16) plt.xlabel("Logarithm of CO2 amount (log10 pCO2 ppmvol)", fontsize=12) plt.ylabel("Mean temperature degC", fontsize=12) plt.axhline(y=100, linestyle="--", color="blue", lw=2, label="Water boils") plt.axhline(y=0, linestyle="--", color="blue", lw=2, label="Water freezes")

plt.scatter(math.log10(280), 13.8, s=200, marker="o", color="green")

plt.xticks(fontsize=12) plt.yticks(fontsize=12)

plt.show()

quit(-1)

1. 1. plotting routines

if(plotvar==0):

       num_steps_per_year = int(ebm.time['num_steps_per_year'])
       Tyear = np.empty((ebm.lat.size, num_steps_per_year))
       for m in range(num_steps_per_year):
           ebm.step_forward()
           Tyear[:,m] = np.squeeze(ebm.Ts)
       Tmin=round(np.min(Tyear),1)
       Tmax=round(np.max(Tyear),1)
       
       #Tmean=round( np.mean(Tyear),1)
       tmeans1=np.mean(Tyear, axis=1)
       #print (ebm.lat)
       latrads1=np.radians(ebm.lat)
       latcoeffs1=np.power(np.cos(latrads1),2)
       tmeans2=tmeans1*latcoeffs1
       Tmean=np.mean(tmeans2)
       #print (np.shape(tmeans1))
       #quit(-1)        
       fig = plt.figure(figsize=(5,5))
       ax = fig.add_subplot(111)
       
       factor = 365. / num_steps_per_year
       cmap1=plt.cm.seismic
       cmap1=plt.cm.winter
       cmap1=plt.cm.coolwarm
       #cmap1=plt.cm.cool_r
       #cmap1=plt.cm.cool
       #cmap1=cmap1.reversed()      
       #levels1=[-80,-70,-60,-50,-40,-30]
       levels2=[-200,-150,-120,-100,-70,-60,-50,-40,-30,-20,-10,0,5,10,15,20,25,30,35,40,45,50,60,70,80,90,100,120,150,300]
       Tyear2 = resize(Tyear, (Tyear.shape[0] *3, Tyear.shape[1]*2),anti_aliasing=True)
       ax.imshow(Tyear, origin="lower", extent=[0,360,-90,90], cmap=cmap1, interpolation="bicubic")
       #cax = ax.contourf(factor * np.arange(num_steps_per_year),
       #              ebm.lat, Tyear[:,:], 
       #              cmap=cmap1, vmin=Tmin, vmax=Tmax, antialiased=False, levels=levels2)
       cs1 = ax.contour(factor * np.arange(num_steps_per_year), ebm.lat,Tyear[:,:],
                     origin="lower", extent=[0,360,-90,90], alpha=0.5, 
                     colors='#00005f', vmin=Tmin, vmax=Tmax, levels=levels2)
       ax.clabel(cs1, cs1.levels, inline=True, fontsize=14)                     
       #cbar1 = plt.colorbar(cax)
       title1='Temperatures degC of planet, if ecc='+str(round(ecc,3))++str(round(long_peri,2))+' and obliquity='+str(round(obliquity,2))+' deg \n if S0= '+ str(round(S1,1)) +' W m-2 , pressure of CO2= '+str(round(co2,2))+' ppm volume'
       title2="\nTemperature deg C:  min "+str(round(Tmin,2))+" max "+str(round(Tmax,2))+" mean "+str(round(Tmean,2))
       #ax.set_suptitle(title1, fontsize=12)
       ax.set_title(title1+title2, fontsize=11)
       ax.tick_params(axis='x', labelsize=12)
       ax.tick_params(axis='y', labelsize=12)
       ax.set_xlabel('Days of year', fontsize=13)
       ax.set_ylabel('Latitude', fontsize=13)
       plt.savefig('1000dpi.svg', dpi=1000)

if(plotvar==1):

       if 'Tf' in ebm.subprocess['albedo'].param.keys():
           Tf = ebm.subprocess['albedo'].param['Tf']
       else:
           print('No ice considered in this model. Can not plot.')

       num_steps_per_year = int(ebm.time['num_steps_per_year'])
       ice_year = np.empty((ebm.lat.size, num_steps_per_year))
       for m in range(num_steps_per_year):
           ebm.step_forward()
           ice_year[:,m] = np.where(np.squeeze(ebm.Ts) <= Tf, 0, 1)
       
       fig = plt.figure(figsize=(5,5))
       ax = fig.add_subplot(111)
       
       factor = 365. / num_steps_per_year
       cax = ax.contourf(factor * np.arange(num_steps_per_year),
                     ebm.lat, ice_year[:,:], 
                     cmap=plt.cm.seismic, vmin=0, vmax=1, levels=2)
       cbar1 = plt.colorbar(cax)
       
       ax.set_title('Ice throughout the year', fontsize=14)
       ax.set_xlabel('Days of year', fontsize=14)
       ax.set_ylabel('Latitude', fontsize=14)

if(plotvar==2):

       fig = plt.figure(figsize=(7.5,5))

       # Temperature plot
       ax2 = fig.add_subplot(111)
       ax2.plot(ebm.lat,ebm.albedo)

       ax2.set_title('Albedo', fontsize=14)
       ax2.set_xlabel('latitude', fontsize=10)
       ax2.set_ylabel(, fontsize=12)

       ax2.set_xticks([-90,-60,-30,0,30,60,90])
       ax2.set_xlim([-90,90])
       ax2.set_ylim([0,1])
       ax2.grid()
   
       plt.show()

plt.suptitle("Planet that has 90 obliquity") plt.title("Temperature deg C") plt.legend() plt.show()