1. 1. temperatures, if S=0.9*S0
   2. sun radiation down 10% from current.
2. python3/climblab code
4. 10.5.2023 0000.0001

import numpy as np import matplotlib.pyplot as plt from matplotlib import cm import climlab from climlab import constants as const

def plot_temp_section(model, timeave=True):

   fig = plt.figure()
   ax = fig.add_subplot(111)
   #viridis = cm.get_cmap('jet')
   #viridis = cm.get_cmap('turbo')
   #viridis = cm.get_cmap('winter')
   viridis = cm.get_cmap('cool_r')
   #viridis = cm.get_cmap('PuBu')
   plt.set_cmap(viridis)
   if timeave:
       field = model.timeave['Tatm'].transpose()
   else:
       field = model.Tatm.transpose()
   cax = ax.contourf(model.lat, model.lev,field-273.15)
   CS = ax.contour(model.lat, model.lev,field-273.15,
                colors='k',  # negative contours will be dashed by default
                )
   ax.clabel(CS,fmt='%1.1f',fontsize=14, inline=1)
   ax.invert_yaxis()
   ax.set_title("Temperature profile, if S=0.93*S0", fontsize=18)
   ax.set_xlabel("Latitude", fontsize=15)
   ax.set_ylabel("Pressure", fontsize=15)
   ax.xaxis.set_tick_params(labelsize=14)
   ax.yaxis.set_tick_params(labelsize=14)
   ax.set_xlim(-90,90)
   ax.set_xticks([-90, -60, -30, 0, 30, 60, 90])
   cbar1=fig.colorbar(cax)
   cbar1.ax.tick_params(labelsize=15)

rau=1.0 ## planet a au S1=1365.2

1. insok=1/(rau*rau) ## insolation coefficient"

insok=0.93 alb=0.299

1. alb=0.06

greenhouse=0.0 cloudiness=1 waterdepth=100

print(rau, insok)

1. 1. not used

delta_t = 60. * 60. * 24. * 30

absorber_vmr = {'CO2':420/1e6,

               'CH4':0.,
               'N2O':0.,
               'O2':0.,
               'CFC11':0.,
               'CFC12':0.,
               'CFC22':0.,
               'CCL4':0.,
               'O3':0.}

1. state = climlab.column_state(num_lev=20, num_lat=1, water_depth=5.)

state = climlab.column_state(num_lev=12, num_lat=16, water_depth=waterdepth)

insol = climlab.radiation.DailyInsolation(name='Insolation',

                                         domains=state['Ts'].domain, S0=S1*insok)

1. olr = climlab.radiation.Boltzmann(name='OutgoingLongwave',state=state, tau = 0.612,eps = 1.,timestep = delta_t)

1. asr = climlab.radiation.SimpleAbsorbedShortwave(name='AbsorbedShortwave',
2. state=state,
3. insolation=341.3,
4. insolation=insol.insolation,
5. albedo=alb,
6. timestep = delta_t)

1. rcm=climlab.TimeDependentProcess(state=state)

1. h2o = climlab.radiation.ManabeWaterVapor(name='H2O', state=state)

h2o = climlab.radiation.ManabeWaterVapor(state=state, relative_humidity=0.5)

1. CAM3 radiation with default parameters and interactive water vapor

1. rad = climlab.radiation.CAM3(name='Radiation', state=state,specific_humidity=h2o.q, albedo=alb)

rad = climlab.radiation.CAM3(name='Radiation', state=state, return_spectral_olr=True, icld=cloudiness, S0 = insol.S0*insok*(1+greenhouse), insolation=insol.insolation, coszen=insol.coszen

1. absorber_vmr = absorber_vmr2

)

print(insol.S0)

1. rad = climlab.radiation.CAM3(name='Radiation',
2. state=state,
3. specific_humidity=h2o.q,
4. S0 = insol.S0,
5. insolation=insol.insolation,
6. coszen=insol.coszen)

1. rad = climlab.radiation.RRTMG_LW(state=state,
2. specific_humidity=h2o.q,
3. S0 = insol.S0,
4. insolation=insol.insolation,
5. icld=0, # Clear-sky only!
6. return_spectral_olr=False, # Just return total OLR
7. absorber_vmr = absorber_vmr)

conv = climlab.convection.ConvectiveAdjustment(name='Convective Adjustment',state=state, adj_lapse_rate=6.5)

rcm = climlab.couple([rad,conv,h2o, insol], name='RCM')

1. print(rcm)

1. quit(-1)

1. rcm.add_subprocess('Radiation', rad)
2. rcm.add_subprocess('WaterVapor', h2o)
3. rcm.add_subprocess('Convection', conv)

1. rcm.integrate_years(1)

rcm.integrate_years(1)

1. fig, ax = plt.subplots(dpi=100)

1. state['Tatm'].to_xarray().plot(ax=ax, y='lev', yincrease=False)

1. state['Tatm'].to_xarray().plot(ax=ax,x='lat', y='lev', yincrease=False)

tatm=state['Tatm']-273.15

1. quit(-1)

1. Create and exact clone of the previous model

diffmodel = climlab.process_like(rcm)

diffmodel.name = 'Seasonal RCE with heat transport'

1. thermal diffusivity in W/m**2/degC

D = 0.05

1. D=0.0001

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

1. diffmodel = climlab.couple([rad,conv,h2o, insol,d], name='Seasonal diffmodel')

print(diffmodel)

diffmodel.integrate_years(1) diffmodel.integrate_years(50)

tatm2=state['Tatm']-273.15

print(tatm2)

print("Plot ")

plot_temp_section(rcm, timeave=True)

1. plot_temp_section(diffmodel, timeave=True)

plot_temp_section(diffmodel, timeave=True)

tlayer1=tatm[...,11].ravel() tlayer2=tatm[...,11].ravel()

1. print (" Tatmlen",len(tlayer1))

tlayer1=np.nan_to_num(tlayer1) tlayer2=np.nan_to_num(tlayer2)

meantemp=np.mean(tlayer1) meantemp2=np.mean(tlayer2)

print(tlayer1) print(tlayer2)

print(" meantemp A ",meantemp) print(" meantemp B ",meantemp2)

1. plot_temp_section(rcm, timeave=True)

1. plt.imshow(tatm)

1. ax.set_xlabel("Temperature (K)")
2. ax.set_ylabel("Pressure (hPa)")
3. ax.grid()

1. plt.plot()

plt.show()