File:Earthlike ocean planet mean temperature difference as function of obliquity and eccentricity 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/climlab code
5. 15.11.2023 0000.0005a2

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=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) 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) return(Tdelta_year)

Sk=1.0 albedo=0.28 co2=280 ecc=0.0 long_peri=0

1. obliquity=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]

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

Tss0=[]

lenum=len(eccentricities0) lenun=len(obliquities0)

for m in range(0,lenum):

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

eccentricities=np.array(eccentricities0) obliquities=np.array(obliquities0) Tst1=np.array(Tss0) Tst=Tst1.reshape(lenun, lenum)

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

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

plt.imshow(Tst2, origin="lower",cmap="coolwarm", vmin=0, vmax=300)

1. plt.imshow(Tst2, origin="lower",cmap="coolwarm", vmin=0, vmax=50)
2. plt.imshow(Tst, origin="lower",cmap="coolwarm", interpolation="bicubic", vmin=0, vmax=50)

cs=plt.contour(Tst2, origin="lower", levels=[-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. csf=plt.contourf(Tst2, origin="lower", levels=[0,10,15,20,25,30,35,40,45,50], colors=["blue", "green","yellow","orange", "red"], alpha=0.5)
2. csf=plt.contourf(Tst2, origin="lower", levels=[0,10,15,20,25,30,35,40,45,50], cmap="coolwarm")

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

1. plt.yticks([0,90])

1. plt.yticks([0,90])

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

1. cs.labels()

1. plt.plot(obliquities, Tst, lw=3, color="#7f0000")
2. plt.title("Earth-like ocean planet EBM \n Obliquity ... T_mean degC", fontsize=16)
3. plt.xlabel("Obliquity", fontsize=14)
4. plt.ylabel("Mean temperature degC", fontsize=14)
5. plt.axhline(y=100, linestyle="--", color="blue", lw=2, label="Water boils")

1. plt.axhline(y=0, linestyle="--", color="blue", lw=2, label="Water freezes")

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

1. plt.yticks(eccentricities, fontsize=14)
2. plt.xticks(obliquities, fontsize=14)

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 Temperature global, year max difference Tmax-Tmin \n as function of obliquity and eccentricity")

ax.set_yticks([0,10,20,30,40,50,60,70,80,90])

ax.set_yticklabels(['0.0', '0.1','0.2', '0.3', '0.4','0.5', '0.6', '0.7','0.8','0.9'], fontsize=12) ax.set_xticks([0,10,20,30,40,50,60,70,80,90]) ax.set_xticklabels(['0', '10','20', '30', '40','50', '60', '70','80','90'], fontsize=12)

ax.set_ylabel("Eccentricity e", fontsize=12) ax.set_xlabel("Obliquity degrees", fontsize=12)

plt.show()

quit(-1)