2. 1. snowball or slushball earth
3. 1. climlab python3 source code
5. 24.10.2023 0000.0003c1

import numpy as np import matplotlib.pyplot as plt import xarray as xr import climlab from climlab import constants as const from climlab.dynamics import MeridionalDiffusion

numyears=50

S1_now=1361.5 ## current sol

1. S1=1*0.93
2. S1=1*0.416 ## Sucellus planet ts Early mars-like
3. S1=1357.18*0.7 ## cretaceous impact
   1. 0.93 snowy ball
4. S1=1365.2*0.82 ## huronian
5. S1=0.77
6. S1=1/1
7. S1=0.44
8. S1=1/1.15

1. S1=0.93
2. S1=0.83 ## huronian

1. https://www.researchgate.net/publication/2225553_Solar_Models_Current_Epoch_and_Time_Dependences_Neutrinos_and_Helioseismological_Properties

S1=0.939 ## 3.8 Ga from oringin S1=0.954 ## 4.0 Ga S1=0.844 ## 2.4 Ga S1=0.677 # 0 Ga S1=0.721 ## 0.2 Ga

1. S1=0.744 ##0.6 ga

S1=0.7

albedo=0.25

waterdepth1=5 cloudiness=0.5

S1_abs=S1_now*S1

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

1. 1. faint young sun ch4++

1. 1. co2=280e-6
2. co2=1
3. o3=1/1e6
4. ch4=800/1e9
5. no2=270/1e9
6. o2=1-co2-ch4-no2-o3 ## O2, simulate N2

1. co2=280e-6
2. co2=0.1
3. o3=1/1e6
4. ch4=800/1e9
5. no2=270/1e9
6. o2=1-co2-ch4-no2-o3 ## O2, simulate N2

1. co2=280e-6

co2=0.95 o3=1/1e6 ch4=800/1e9 no2=270/1e9 o2=1-co2-ch4-no2-o3 ## O2, simulate N2

title0=' Earth and faint young sun ' title1='Temperatures throughout the year °C \n if S0 = '+ str(S1) +'*Ssol , CO2 = '+str(round(co2*1e2,2))+' % , CH4 = '+str(round(ch4*1e6,2))+' ppmv'

print(title1)

1. quit(-1)

1. 1. NOTE

absorber_vmr1 = {'CO2':co2,

               'CH4':ch4,
               'N2O':no2,
               'O2':o2,
               'CFC11':1./1e9,
               'CFC12':1./1e9,
               'CFC22':1./1e9,
               'CCL4':1./1e9,
               'O3':o3}

1. 1. ...

num_lev = 24 ## 50

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 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=[-200,-100,-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,200,300,500] Tminv=-100 Tmaxv=120 #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[:,:], 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)

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. main code

state = climlab.column_state(num_lev=num_lev, num_lat=90, water_depth=waterdepth1) lev = state.Tatm.domain.axes['lev'].points

1. 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)

1. indices
2. high = 10 # corresponds to 210 hPa
3. low = 40 # corresponds to 810 hPa

high=1*2 low=9*2

1. A high, thin ice layer (cirrus cloud)
2. r_ice[:,high] = 14. # Cloud ice crystal effective radius (microns)
3. ciwp[:,high] = 10. # in-cloud ice water path (g/m2)
4. cldfrac[:,high] = 0.322
5. A low, thick, water cloud layer (stratus)
6. r_liq[:,low] = 14. # Cloud water drop effective radius (microns)
7. clwp[:,low] = 100. # in-cloud liquid water path (g/m2)
8. cldfrac[:,low] = 0.21

1. 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

1. A low, thick, water cloud layer (stratus)

r_liq[:,low] = 14. # Cloud water drop effective radius (microns) clwp[:,low] = 4. # in-cloud liquid water path (g/m2) cldfrac[:,low] = 0.05

1. wrap everything up in a dictionary

mycloud = {'cldfrac': cldfrac,

         'ciwp': ciwp,
         'clwp': clwp,
         'r_ice': r_ice,
         'r_liq': r_liq}

1. plt.plot(cldfrac[0,:], lev)
2. plt.gca().invert_yaxis()
3. plt.ylabel('Pressure hPa')
4. plt.xlabel('Cloud fraction')
5. plt.title('Cloud fraction in the column model')
6. plt.show()

1. 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)

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)

1. 1. no clouds !!!

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

model.add_subprocess('Radiation', rad) model.add_subprocess('Insolation', sun) model.add_subprocess('WaterVapor', h2o) model.add_subprocess('Convection', conv)

1. print(model.subprocess['Radiation'].state)

1. quit(-1)

1. 1. thermal diffusivity SI units

D = 0.04

1. 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.add_subprocess('SHF', shf) model.add_subprocess('LHF', lhf)

1. 1. model.subprocess['LHF'].Cd *= 0.5

1. 1. 2x co2

1. model.subprocess['LW'].absorptivity = model.subprocess['LW'].absorptivity*1.1

1. quit(-1)

1. One more year to get annual-mean diagnostics
2. model.step_forward()

1. model.integrate_years(1.)

model.integrate_years(numyears)

1. 1. this is very slooow ...
2. model.integrate_converge()

lat = model.lat

num_steps_per_year = int(model.time['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)

plotmonths2(model, Tss-273.15, lat)

plt.show()

print(".")

quit(-1)