File:Ocean planet insolation and co2 effect to average temperature 1 1 1 1.png

https://climlab.readthedocs.io/en/latest/
https://github.com/climlab/climlab

https://paleodyn.uni-bremen.de/study/Climlab_lecture_solution.html

An Introduction to Climlab

Author of this notebook: Maybritt Schillinger

https://climlab.readthedocs.io/en/latest/courseware/Seasonal_cycle_and_heat_capacity.html

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/climlab code
5. 16.11.2023 0000.0006ax1

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 import skimage 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_01(Sk, albedo0, co2, ecc, long_peri, obliquity): numyears=30 ##n no function here, run to equil numlat=18 numlev=6 plotvar=0 ## 1,2,3 lot temp, ice, mean albedo waterdepth=70 S1=1361.5*Sk #print(S1)

diffu1=0.3 # meridional diffusivity in m**2/s #diffu1=1.0 # meridional diffusivity in m**2/s #diffu1=0.01 #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=-100, **ebm.param) #alb = climlab.surface.albedo.ConstantAlbedo(domains=surface, **ebm.param)

#alb.albedo=alb.albedo*0+albedo #quit(-1)

#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.Ts.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) min_year = np.empty(num_steps_per_year) max_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() min_year[m] = np.min(ebm.Ts) max_year[m] = np.max(ebm.Ts)

Tmean_year = np.mean(mean_year) Tmin_year = np.mean(min_year) Tmax_year = np.mean(max_year) Tdelta_year = Tmax_year-Tmin_year #print(round(Tmean_year,2)) #return(Tmean_year) #print(ebm.albedo) #plt.imshow(ebm.Ts) #plt.show() return(Tmean_year)

Sk=1.0

1. albedo=0.28
2. albedo=0.33

albedo=0.15 co2=280 ecc=0.0 long_peri=0

obliquity=0.0

1. eccentricities0=[0,0.9]

1. obliquities0=[0,90]

1. obliquities0=[0,30,60,90]
2. eccentricities0=[0,0.3,0.6,0.9]

1. obliquities0=[0,10,20,30,40,50,60,70,80,90]
2. eccentricities0=[0,.10,.20,.30,.40,.50,.60,.70,.80,.90]

1. Sks0=[0.5,1,2,4]
2. co2s0=[1e-6,280e-6,0.1,1]

1. Sks0=[0.5,1,2,4]
2. co2s0=[1e-6,280e-6,0.1,1]

Sks0=[0.4,0.6, 0.8,1.0,1.2,1.4,1.6] co2s0=[1, 10, 100,1000,1e4, 1e5, 1e6]

1. Sks0=[0.5,2]
2. co2s0=[1e-6,0.1]

1. Ts=run_ebm_01(Sk, albedo, co2, ecc, long_peri, obliquity)

1. print(Sk, co2, Ts)

1. quit(-1)

Tss0=[]

lenum=len(Sks0) lenun=len(co2s0)

for m in range(0,lenum):

   Sk=Sks0[m]
   for n in range(0,lenun):
       co2=co2s0[n]
       Ts=run_ebm_01(Sk, albedo, co2, ecc, long_peri, obliquity)
       print(Sk, co2, Ts)
       Tss0.append(Ts)

Sks=np.array(Sks0) co2s=np.array(co2s0) Tst1=np.array(Tss0) Tst=Tst1.reshape(lenun, lenum)

1. 1. Tst2 = skimage.transform.resize(Tst, (70, 70), anti_aliasing=True, anti_aliasing_sigma=1)

Tst2 = skimage.transform.resize(Tst, (70, 70), anti_aliasing=True)

1. plt.imshow(Tst, origin="lower",cmap="coolwarm", vmin=-200, vmax=200)

plt.imshow(Tst2, origin="lower",cmap="coolwarm", interpolation="bicubic", vmin=-200, vmax=200)

cs=plt.contour(Tst2, origin="lower", levels=[-250,-200,-150,-120,-100,-80,-50,-25,-10,0,10,15,20,25,30,35,40,45,50,60,70,80,90,100,150,200,300,400,500,600,700,800,1000], colors=["#00003f"], alpha=0.5)

1. cs=plt.contour(Tst, origin="lower", levels=[0,10,15,20,25,30,35,40,45,50,75,100], color="#3f0000", alpha=0.5)

plt.clabel(cs, cs.levels, inline=True, fmt=f"%.1f", fontsize=10)

1. cs.labels()

ax=plt.gca()

1. labels1=[item.get_text() for item in axes.get_xticklabels()]

1. xlabels1=["0.",""]

1. axes.set_xticklabels(xlabels1)

plt.title("Ocean planet \n Global average temperature degC \n as function of insolation and pCO2")

p1x=30 p1y=int(math.log10(280)*10)

ax.scatter(p1x,p1y, s=200, marker="o", color="green", alpha=0.3)

1. ax.set_yticks([0,1,2,3])
2. ax.set_yticklabels(['1ppm', '280ppm',"0.1", "1"], fontsize=12)

ax.set_yticks([0,10,20,30,40,50,60]) ax.set_yticklabels(["1 ppm", "10 ppm", "100 ppm","1000 ppm","1%", "10%", "100%"])

1. ax.set_xticks([0,1,2,3])

ax.set_xticks([0,10,20,30,40,50,60]) ax.set_xticklabels(["0.4","0.6", "0.8","1.0","1.2","1.4","1.6"])

1. ax.set_xticklabels(['0.5', '1.0','2.0',"4.0"], fontsize=12)

ax.set_xlabel("Insolation S/Sol", fontsize=12) ax.set_ylabel("pCO2 ppmv", fontsize=12)

for n in range(len(Sks)): for m in range(len(co2s)): print(round(Sks[n],3),round(co2s[m]),round(Tst1[m*len(co2s)+n],2))

plt.show()

print(".")

quit(-1)