File:Ocean planet distance vs temperature moist energy balance model 1 1 1 1.png

1. must install additional packages along climlab
3. pip install climlab
4. 1. and source pip install "pip install ." in source directory
6. https://github.com/climlab/climlab-cam3-radiation
7. https://github.com/climlab/climlab-emanuel-convection
8. https://github.com/climlab/climlab-rrtmg

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. orbital distance vs temperature of ocean planet
   1. seasonal climlab energy balance model
3. python3/climlab code
5. must install additional packages along climlab
7. pip install climlab
8. 1. and source pip install "pip install ." in source directory
10. https://github.com/climlab/climlab-cam3-radiation
11. https://github.com/climlab/climlab-emanuel-convection
12. https://github.com/climlab/climlab-rrtmg
14. 17.11.2023 0000.0006a1

import math import numpy as np import matplotlib.pyplot as plt from matplotlib import cm

from scipy.interpolate import griddata import skimage from skimage.transform import resize

import xarray as xr

import climlab from climlab import constants as const from climlab.process.diagnostic import DiagnosticProcess from climlab.domain.field import Field, global_mean from climlab.dynamics import MeridionalDiffusion

def plotmonths(Ts, lat): global title1 lela=len(lat) print(np.shape(Ts)) fig = plt.figure( figsize=(8,5) ) ax = fig.add_subplot(111) clevels=10 Tmin=-50 Tmax=50 plt.xticks(fontsize=15) plt.yticks(fontsize=15) cax = ax.contourf(np.arange(365)+0.5, lat, Ts,cmap=plt.cm.coolwarm,vmin=Tmin, vmax=Tmax, levels=256 ) cc = ax.contour(np.arange(365)+0.5, lat, Ts, colors=['#00003f'],) ax.clabel(cc, cc.levels, colors=['#00005f'], inline=True, fmt='%3.1f',fontsize=15) #ax.set_tick_params(axis='both', which='minor', labelsize=15) ax.set_xlabel('Day', fontsize=15) ax.set_ylabel('Latitude', fontsize=15) #cbar = plt.colorbar(cax) #cbar.set_clim(-50.0, 50.0) ax.set_title('Zonal mean surface temperatures (degC)', fontsize=16) return(0)

def plotmonths2(model, Ts, lat): global title1

#num_steps_per_year=360

Tmin=np.min(Ts) Tmax=np.max(Ts) fig = plt.figure(figsize=(5,5)) ax = fig.add_subplot(111) factor = 365. / num_steps_per_year #cmap1=plt.cm.seismic #cmap1=plt.cm.turbo cmap1=plt.cm.coolwarm #cmap1=plt.cm.winter #cmap1=plt.cm.cool_r #cmap1=plt.cm.cool #cmap1=cmap1.reversed() #levels1=[-80,-70,-60,-50,-40,-30] levels2=[-250,-200,-150,-120,-100,-90,-80,-70,-60,-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50,55,60,65,80,100,120,140,160,180,200,300,500] Tminv=-150 Tmaxv=150 #cax = ax.contourf(factor * np.arange(num_steps_per_year), # ebm.lat, Tyear[:,:], # cmap=cmap1, vmin=Tminv, vmax=Tmaxv, antialiased=False, levels=256) ax.imshow(Ts[:,:],origin="lower", extent=[0,360,-90,90],cmap=cmap1, vmin=Tminv, vmax=Tmaxv, interpolation="bicubic") cs1 = ax.contour(factor * np.arange(num_steps_per_year),model.lat, Ts[:,:], alpha=0.5, origin="lower", extent=[0,360,-90,90],colors='#00005f', vmin=Tminv, vmax=Tmaxv, levels=levels2) ax.clabel(cs1, cs1.levels, inline=True, fontsize=14) #cbar1 = plt.colorbar(cax) ax.set_title(title1, fontsize=12) fig.suptitle(title0, fontsize=22) ##ax_set_suptitle(title0, fontsize=18) 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.tight_layout() plt.savefig('1000dpi.svg', dpi=1000) return(0)

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 model_ebm_00(Sk,ecc, long_peri, obliquity,albedo0,co20,waterdepth,cloudiness): num_years=30 num_lev=4 num_lat=36 num_steps_per_year=72

plotvar=0 ## 1,2,3 lot temp, ice, mean albedo #waterdepth=70 S1=1361.5*Sk

co2=co20 #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=num_lat, lum_lon=None, num_lev=num_lev,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)

#Tf=-10 Tf=-10 alb = climlab.surface.albedo.StepFunctionAlbedo(state=ebm.state, Tf=Tf, **ebm.param) #alb = climlab.surface.albedo.StepFunctionAlbedo(state=ebm.state, Tf=Tf, **ebm.param) #alb = climlab.surface.albedo.ConstantAlbedo(domains=surface, **ebm.param) #albedo=tanalbedo(ebm.state, **ebm.param)

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

#alb = tanalbedo(state=ebm.state, **ebm.param) ebm.remove_subprocess('albedo') 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

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.remove_subprocess('diffusion') ebm.add_subprocess('diffusion', diff) #print (ebm) ebm.time['num_steps_per_year']=num_steps_per_year ebm.step_forward() #ebm.diagnostics #ebm.integrate_years(numyears) ebm.integrate_years(num_years) ##ebm.integrate_converge() #print(ebm.Ts) #plt.plot(ebm.Ts) #plt.show() ebm.time['num_steps_per_year']=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

#temps_greater_0 = (data > 10).sum()

#print(round(Tmean_year,2)) #return(Tmean_year) #print(ebm.albedo) #plt.imshow(ebm.Ts) #plt.show() return(Tmean_year, Tmin_year ,Tmax_year )

1. 1. nok ??

def model_cam3_00(insok,ecc, long_peri, obliquity,albedo,pco2,waterdepth,cloudiness): ############## main code ## note cloudiness not func num_years=10 num_lev=8 num_lat=18 num_steps_per_year=36

S1_base=1361.5 ## current # thermal diffusivity in W/m**2/degC #D = 0.05 relative_humidity=1.00*0.001 #D=0.0001 D=0.5 co2=pco2 S1_abs=S1_base*insok orbit1={'ecc': ecc, 'long_peri': long_peri, 'obliquity': obliquity, 'S0':S1_abs} #delta_t = 60. * 60. * 24. * 60

absorber_vmr = {'CO2':pco2, 'CH4':800./1e9, 'N2O':100./1e9, 'O2':0.21, 'CFC11':1./1e9, 'CFC12':1./1e9, 'CFC22':1./1e9, 'CCL4':1./1e9, 'O3':1./1e6}

#state = climlab.column_state(num_lev=20, num_lat=1, water_depth=5.) state = climlab.column_state(num_lev=num_lev, num_lat=num_lat, water_depth=waterdepth) insol = climlab.radiation.DailyInsolation(name='Insolation', domains=state['Ts'].domain, S0=S1_abs, orb=orbit1) #insol.S0=S1_abs h2o = climlab.radiation.ManabeWaterVapor(state=state, relative_humidity=relative_humidity) rad = climlab.radiation.CAM3(name='Radiation', state=state, return_spectral_olr=True, icld=cloudiness, S0 = S1_abs, insolation=insol.insolation, coszen=insol.coszen, absorber_vmr = absorber_vmr,albedo=albedo) #print(insol.S0) #print (insol.coszen) #quit(-1) conv = climlab.convection.ConvectiveAdjustment(name='Convective Adjustment',state=state, adj_lapse_rate=6.5) #rcm = climlab.couple([insol, rad,conv,h2o], name='RCM') rcm = climlab.couple([insol, rad], name='RCM') surface = rcm.domains['Ts'] rcm.water_depth=waterdepth #rcm.Tf=-10 rcm.Tf=-10 #quit(-1) #rad.a0=albedo ## WARNING DYNAMIC ALBEDO NOK #rcm.remove_subprocess('albedo') #alb = climlab.surface.albedo.StepFunctionAlbedo(state=rcm.state, Tf=-10, **rcm.param) alb = climlab.surface.albedo.ConstantAlbedo(domains=surface, **rcm.param) #alb.albedo[:]=albedo #alb = tanalbedo(state=rcm.state, **rcm.param) rcm.add_subprocess('albedo', alb) #print (alb.diagnostics) #quit(-1) #print(" Integrate ...") rcm.time['num_steps_per_year']=num_steps_per_year

rcm.integrate_years(3) # Create and exact clone of the previous model diffmodel = climlab.process_like(rcm) diffmodel.name = 'Seasonal RCE with heat transport' # meridional diffusivity in m**2/s K = D / diffmodel.Tatm.domain.heat_capacity[0] * const.a**2 #print("K ", K) d = climlab.dynamics.MeridionalDiffusion(K=K, state={'Tatm': diffmodel.Tatm}, **diffmodel.param) #diffmodel.add_subprocess('Meridional Diffusion', d) #diffmodel = climlab.couple([rad,conv,h2o, insol,d], name='Seasonal diffmodel') #print(diffmodel) #diffmodel.integrate_years(1) diffmodel.integrate_years(num_years) #diffmodel.integrate_converge() tatm2=state['Tatm']-273.15 #print(tatm2) #print("Plot ") #plot_temp_section(rcm, timeave=True) #plot_temp_section(diffmodel, timeave=True) #plot_temp_section(diffmodel, timeave=True) tlayer1=tatm2[...,7].ravel() tlayer2=tatm2[...,0].ravel() #print (" Tatmlen",len(tlayer1)) tlayer1=np.nan_to_num(tlayer1) tlayer2=np.nan_to_num(tlayer2) ## check this meantemp=np.mean(tlayer1) meantemp2=np.mean(tlayer2) #print(tlayer1) #print(tlayer2) #print(" meantemp A ",meantemp) #print(" meantemp B ",meantemp2) #fig, ax = plt.subplots(dpi=100) #state['Tatm'].to_xarray().plot(ax=ax, y='lev', yincrease=False) #state['Tatm'].to_xarray().plot(ax=ax,x='lat', y='lev', yincrease=False) tatm=state['Tatm']-273.15 rcm=diffmodel years=1 num_steps_per_year = int(rcm.time['num_steps_per_year']) Tyear = np.empty((rcm.lat.size, num_steps_per_year*years)) Tmean1=0 #print (rcm.lat) cosu=np.cos(np.radians(rcm.lat)) for m in range(num_steps_per_year*years): rcm.step_forward() #Tyear[:,m] = np.squeeze(rcm.Ts) #Tyear[:,m] = np.squeeze(rcm.Ts) tatm2=state['Tatm']-273.15 tlayer1=tatm2[...,7] tlayer1=np.nan_to_num(tlayer1) Tyear[:,m] = np.squeeze(tlayer1) Tmean0=np.mean((rcm.Ts-273.15)*cosu) #Tmean0=np.mean((rcm.Ts-273.15)) #Tmean0=np.mean(tlayer1) Tmean1=Tmean1+Tmean0 #print(cosu) #print(rcm.Ts) #print(Tmean0) #quit(-1)

#Tmean=Tmean1/(num_steps_per_year*years) Tmean=np.mean(Tyear) Tmin=np.min(Tyear) Tmax=np.max(Tyear)

return(Tmean, Tmin, Tmax)

1. 1. nok ?

def model_rrtmg_00(Sok,ecc, long_peri, obliquity,albedo,pco2,waterdepth,cloudiness):

numyears=5

num_lev = 8 ## 50 num_lat=18 num_steps_per_year=36

S1_now=1361.5 ## current sol co2=pco2*1e-6 o3=1/1e6 ch4=800/1e9 no2=270/1e9 o2=1-co2-ch4-no2-o3 ## O2, simulate N2

S1=Sok S1_abs=S1_now*S1

orbit1={'ecc': ecc, 'long_peri': long_peri, 'obliquity': obliquity, 'S0':S1_abs}

absorber_vmr1 = {'CO2':co2, 'CH4':ch4, 'NO2':no2, 'O2':o2, 'CFC11':1./1e9, 'CFC12':1./1e9, 'CFC22':1./1e9, 'CCL4':1./1e9, 'O3':o3}

state = climlab.column_state(num_lev=num_lev, num_lat=num_lat, water_depth=waterdepth) lev = state.Tatm.domain.axes['lev'].points # Define two types of cloud, high and low cldfrac = np.zeros_like(state.Tatm) r_liq = np.zeros_like(state.Tatm) r_ice = np.zeros_like(state.Tatm) clwp = np.zeros_like(state.Tatm) ciwp = np.zeros_like(state.Tatm) # indices #high = 10 # corresponds to 210 hPa #low = 40 # corresponds to 810 hPa high=1 low=7 # A high, thin ice layer (cirrus cloud) #r_ice[:,high] = 14. # Cloud ice crystal effective radius (microns) #ciwp[:,high] = 10. # in-cloud ice water path (g/m2) #cldfrac[:,high] = 0.322 # A low, thick, water cloud layer (stratus) #r_liq[:,low] = 14. # Cloud water drop effective radius (microns) #clwp[:,low] = 100. # in-cloud liquid water path (g/m2) #cldfrac[:,low] = 0.21 # A high, thin ice layer (cirrus cloud) r_ice[:,high] = 14. # Cloud ice crystal effective radius (microns) ciwp[:,high] = 2. # in-cloud ice water path (g/m2) cldfrac[:,high] = 0.1 # A low, thick, water cloud layer (stratu r_liq[:,low] = 14. # Cloud water drop effective radius (microns) clwp[:,low] = 4. # in-cloud liquid water path (g/m2) cldfrac[:,low] = 0.05 # wrap everything up in a dictionary mycloud = {'cldfrac': cldfrac, 'ciwp': ciwp, 'clwp': clwp, 'r_ice': r_ice, 'r_liq': r_liq} #plt.plot(cldfrac[0,:], lev) #plt.gca().invert_yaxis() #plt.ylabel('Pressure hPa') #plt.xlabel('Cloud fraction') #plt.title('Cloud fraction in the column model') #plt.show() #quit(-1) model = climlab.TimeDependentProcess(state=state, name='Radiative-Convective-Diffusive Model') h2o = climlab.radiation.ManabeWaterVapor(state=state) conv = climlab.convection.ConvectiveAdjustment(state={'Tatm':model.state['Tatm']},adj_lapse_rate=6.5,**model.param) sun = climlab.radiation.DailyInsolation(name='Insolation', domains=state['Ts'].domain, S0=S1_abs, orb=orbit1) #print (sun.S0) #sun.S0=S1 #sun.ecc=ecc #sun.long_peri=long_peri #sun.obliquity=obliquity #print (sun.insolation) #quit(-1) #rad = climlab.radiation.RRTMG(state=state, #specific_humidity=h2o.q, #albedo=albedo, #S0=S1_abs, #co2vmr=co2,ch4vmr=ch4, #n2ovmr=no2,o2vmr=o2, #cfc11vmr=0.0,cfc12vmr=0.00,cfc22vmr=0.00, #ccl4vmr=0.0,o3vmr=o3, #insolation=sun.insolation, #coszen=sun.coszen, #absorber_vmr = absorber_vmr1, #**mycloud) #quit(-1) rad = climlab.radiation.RRTMG(state=state,specific_humidity=h2o.q, albedo=albedo, S0=S1_abs, co2vmr=co2,ch4vmr=ch4, n2ovmr=no2, o2vmr=o2, cfc11vmr=0.0,cfc12vmr=0.00,cfc22vmr=0.00, ccl4vmr=0.0,o3vmr=o3, insolation=sun.insolation, orbit=orbit1, coszen=sun.coszen, #absorber_vmr = absorber_vmr1, )

## no clouds !!! #rad = climlab.radiation.CAM3(name='Radiation', state=state,return_spectral_olr=True,icld=cloudiness,S0 = S1_abs, #insolation=sun.insolation,coszen=sun.coszen,albedo=albedo,absorber_vmr = absorber_vmr3)

model.add_subprocess('Radiation', rad) #model.remove_subprocess('Insolation') model.add_subprocess('Insolation', sun) model.add_subprocess('WaterVapor', h2o) model.add_subprocess('Convection', conv) #print(model.subprocess['Radiation'].state) #quit(-1) ## thermal diffusivity SI units D = 0.04 # meridional diffusivity SI units K = D / model.Tatm.domain.heat_capacity[0] * const.a**2 d = MeridionalDiffusion(state={'Tatm': model.state['Tatm']}, K=K, **model.param) #model.add_subprocess('Diffusion', d) shf = climlab.surface.SensibleHeatFlux(state=model.state, Cd=0.5E-3) lhf = climlab.surface.LatentHeatFlux(state=model.state, Cd=0.5E-3) lhf.q = h2o.q #model.remove_subprocess('SHF') #model.remove_subprocess('LHF') model.add_subprocess('SHF', shf) model.add_subprocess('LHF', lhf) #model.subprocess['LHF'].Cd *= 0.5 ## 2x co2 #model.subprocess['LW'].absorptivity = model.subprocess['LW'].absorptivity*1.1 #quit(-1) # One more year to get annual-mean diagnostics #model.step_forward() #model.integrate_years(1.) #quit(-1) model.integrate_years(numyears) #quit(-1) ## this is very slooow ... #model.integrate_converge() lat = model.lat #num_steps_per_year = int(model.time['num_steps_per_year']) model.time['num_steps_per_year']=num_steps_per_year

Tss = np.empty((lat.size, num_steps_per_year)) for n in range(num_steps_per_year): model.step_forward() Ts=model.Ts Tss[:,n] = np.squeeze(Ts)

coslat=np.cos(np.radians(model.lat)) sheip1=np.shape(Tss) meanstime=np.mean(Tss-273.15, axis=1) #print(meanstime)

Tmean=np.mean(coslat*meanstime) Tmin=np.min(Tss-273.15) ## all year mon Tmax=np.max(Tss-273.15) ## all year min

#plotmonths2(model, Tss-273.15, lat) #plt.show() ## WARNING nok, cos lat not accounted! #Tmean=np.mean(Tss-273.15) return(Tmean, Tmin, Tmax)

1. 2. main program

Sk=1.0

1. albedo=0.28
2. albedo=0.33

albedo=0.33

ecc= 0.0167643 long_peri=280.32687 obliquity=23.459277

co2=280

1. ecc=0.0
2. long_peri=0
3. obliquity=0

waterdepth=50

cloudiness=0.0

emiss=0.95 sb=5.670374419e-8

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]

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

1. Sks0=[0.3, 0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1]
2. co2s0=[1, 10, 100,1000,1e4,1e5,1e6,10e6,100e6]

1. Sks0=np.arange(0.2,2.0,0.2).tolist()

aus0=np.arange(0.7,1.3,0.01).tolist()

1. Sks0=[0.5,2]
2. co2s0=[10,0.1*1e6]

1. Sks0=[0.6,1.0, 1.4]
2. co2s0=[10,280,100000]

1. 1. base seasonal a+bt ebm

1. 1. cam3 radiation model

1. 1. rrtmg model

Sk=1.0

1. Sk=0.35

co2=280

1. co2=1e6*100

Ts, Tmin, Tmax =model_ebm_00(Sk,ecc, long_peri, obliquity,albedo,co2,waterdepth,cloudiness)

1. Ts=model_cam3_00(Sk,ecc, long_peri, obliquity,albedo,co2,waterdepth,cloudiness) ## nok
2. Ts=model_rrtmg_00(Sk,ecc, long_peri, obliquity,albedo,co2,waterdepth,cloudiness)

print(Sk, co2, Ts)

1. quit(-1)

Tss0=[] Tmins0=[] Tmaxs0=[] Teqs0=[]

1. aus0=[]

Sks0=[]

lenum=len(aus0)

1. lenum=len(Sks0)
2. lenun=len(co2s0)

co2=280

Kelvin=273.15

for m in range(0,lenum): #Sk=Sks0[m] #au=1/math.pow(Sk,2) au=aus0[m] Sk=1/math.pow(au,2) S1_abs=Sk*1361.5

# for n in range(0,lenun): # co2=co2s0[n]

Teq=(math.pow(Sk, 1/4)*288)-273.15 ## teqk 255 without greenhouse Teq_slow = np.power(((S1_abs*(1-albedo)) / (1*sb)), 1/4)-Kelvin Teq_fast = np.power(((S1_abs*(1-albedo)) / (4*sb)), 1/4)-Kelvin ## three distinct model Ts1, Tmin1, Tmax1=model_ebm_00(Sk,ecc, long_peri, obliquity,albedo,co2,waterdepth,cloudiness) #Ts2, Tmin2, Tmax2=model_cam3_00(Sk,ecc, long_peri, obliquity,albedo,co2,waterdepth,cloudiness) #Ts3, Tmin3, Tmax3=model_rrtmg_00(Sk,ecc, long_peri, obliquity,albedo,co2,waterdepth,cloudiness) #Ts=(Ts1+Ts2+Ts3)/3 #Tmin=(Tmin1+Tmin2+Tmin3)/3 #Tmax=(Tmax1+Tmax2+Tmax3)(3 Ts=Ts1 Tmin=Tmin1 Tmax=Tmax1

print(au, Sk, co2, Ts) #aus0.append(au) Tss0.append(Ts) Tmins0.append(Tmin) Tmaxs0.append(Tmax) Sks0.append(Sk) Teqs0.append((Teq))

Sks=np.array(Sks0)

1. co2s=np.array(co2s0)

Tst1=np.array(Tss0) Tmax1=np.array(Tmaxs0) Tmin1=np.array(Tmins0) aus=np.array(aus0) Teqs1=np.array(Teqs0)

plt.plot(aus, Tst1, lw=3, color="green") plt.plot(aus, Tmin1, lw=2, color="blue") plt.plot(aus, Tmax1, lw=2, color="red") plt.plot(aus, Teqs1, lw=1.5, linestyle=":",color="#3f0000", label="Equil. temp. Earth-like")

plt.axhline(y=0, color="#7f7fff", linestyle=":", lw=2, alpha=0.5, label="Water boils") plt.axhline(y=100, color="blue", linestyle=":",lw=2, alpha=0.5,label="Water freezes") plt.axhline(y=15, color="green", linestyle="--",lw=2, alpha=0.5,label="Earth now") plt.axhline(y=50, color="red", linestyle="dashdot",lw=2, alpha=0.5,label="Human habitability limit?")

plt.title(" Ocean planet, Sun-like star \n mean temperature vs distance AU", fontsize=18) plt.xlabel("Distance AU", fontsize=16) plt.ylabel("Mean temperature °C", fontsize=16) plt.xticks(fontsize=14) plt.yticks(fontsize=14)

plt.legend() plt.grid()

print(aus) print(Tst1) print(Teqs1) print(".") plt.show()

quit(-1)

1. 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, (80, 80), anti_aliasing=False)

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

Tst3=plt.imshow(Tst, origin="lower",cmap="coolwarm", interpolation="bicubic", vmin=-50, vmax=50)

cs=plt.contour(Tst3.data, origin="lower", levels=[-250,-200,-150,-120,-100,-80,-70,-60,-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,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=70 p1y=int(math.log10(280)*10)

1. 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)
3. ax.set_yticks([0,10,20,30,40,50,60])
4. ax.set_yticklabels(["1 ppm", "10 ppm", "100 ppm","1000 ppm","1%", "10%", "100%"])

1. ax.set_xticks([0,1,2,3])
2. ax.set_xticks([0,10,20,30,40,50,60])
3. 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)

1. Sks0=[0.3, 0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1]
2. co2s0=[1, 10, 100,1000,1e4,1e5,1e6,10e6,100e6]

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

ax.set_xticks([0,10,20,30,40,50,60,70,80]) ax.set_xticklabels(["0.3","0.4", "0.5","0.6","0.7","0.8","0.9","1.0","1.1" ])

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)