File:Impact winter ebm simulation 2 mean temperature changes 1 1.svg

###############3
##
## impact winter simulation
## python, climlab, energy balance
##
## 17.5.2023 0000.0002
#
######

import numpy as np
import matplotlib.pyplot as plt
import math
import xarray as xr

import climlab
from climlab import constants as const

years1=15

no_dimming=1

##S1=1365.2*1.0 ## present
S1=1357.18*1.0 ## K-Pg limit

CO2=500 ## CO2 500 ppm

minidim=0.999 ## S, if minimun sun dimming
minilimit=720 ##  days , that minimun sun dimming is possible

waterdepth1=19

#timescale1=2
timescale1=0.5 ## dimming recovert exponential time scale, in years

snapshot_year=2

def plot_tann(Tann, lats):
    lat = model1.lat
    plt.plot(lat, Tann)
    plt.xlim(-90,90)
    plt.xlabel('Latitude', fontsize=15)
    plt.ylabel('Temperature (degC)', fontsize=15)
    plt.title('Annual mean temperature', fontsize=15)
    plt.legend( "Ts" , fontsize=15)
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    plt.show()

def plotmonths(Ts, lat):
    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(4*np.arange(90)+0.5, lat, Ts,cmap=plt.cm.coolwarm,vmin=Tmin, vmax=Tmax, levels=256 )
    cc = ax.contour(4*np.arange(90)+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)
    
def plotmonthcurve(Ts, lat, latitude):
    lela=len(lat)
    print(np.shape(Ts))
    fig = plt.figure( figsize=(8,5) )
    ax = fig.add_subplot(111)
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    pp=int((latitude+90)/2)
    lattari=lat[pp]
    print(latitude, lattari)
    Tk=Ts[pp]
    meant1=np.mean(Tk)
    print("Lat ", latitude, "Mean T annual ", meant1)
    plt.plot(4*np.arange(90)+0.5, Tk)
    ax.set_xlabel('Day', fontsize=15)
    ax.set_ylabel('T', fontsize=15)

    
#model1 = climlab.EBM_seasonal()

#model1 = climlab.EBM_seasonal(D=2.5, num_points = 360, a0=0.3, a2=0.078, ai=0.62, water_depth=5)
model1 = climlab.EBM_seasonal(D=1, num_points = 360, a0=0.29, a2=0.078, ai=0.62, water_depth=waterdepth1)

model1.subprocess['insolation'].S0=S1

#LW= climlab.radiation.Boltzmann(eps=0.65, tau=0.95,state=model1.state,**model1.param)

#LW= climlab.radiation.Boltzmann(eps=0.58, #tau=1.0,state=model1.state,**model1.param)

LW = climlab.radiation.aplusbt.AplusBT_CO2(CO2=CO2, state=model1.state, **model1.param)

model1.remove_subprocess('LW')

model1.add_subprocess('LW',LW)

#budyko1 = climlab.dynamics.BudykoTransport(b=3.81,state=model1.state,**model1.param)

#model1.add_subprocess('budyko_transport',budyko1)

# thermal diffusivity in W/m**2/degC
D = 1e-6
#D=0.0001

# meridional diffusivity in m**2/s
#K = D / model1.Ts.domain.heat_capacity[0] * const.a**2
K=6.5e5
print("K ", K)
diff1 = climlab.dynamics.MeridionalDiffusion(K=K, state={'Ts': model1.Ts}, **model1.param)


model1.remove_subprocess('diffusion')

model1.add_subprocess('diffusion', diff1)

model1.diagnostics.keys()
model1.state.keys()

model1.integrate_days(50.)

model1.integrate_years(1, verbose=True)

##model1.icelat
Tequil0=np.array(model1.Ts)
Tequil0_mean = climlab.global_mean(model1.Ts)

Tequil0_min = np.min(Tequil0)
Tequil0_max = np.max(Tequil0)

print(Tequil0_mean,Tequil0_min,Tequil0_max )
print(Tequil0_max-Tequil0_min)

#model1.OLR=model1.OLR*0.5

lat = model1.lat

num_steps_per_year = int(model1.time['num_steps_per_year'])
#print(num_steps_per_year)

## integrate to equilibrium 
model1.integrate_converge()

Tequil = np.array(model1.Ts)
Tequil_mean = climlab.global_mean(model1.Ts)
Tequil_min = np.min(Tequil)
Tequil_max = np.max(Tequil)

print(Tequil_mean,Tequil_min,Tequil_max )

print(Tequil_mean-Tequil_mean)

#quit(-1)

#Tss=[]

Tss = np.empty((lat.size, num_steps_per_year*years1))

tmeans1=[]
Sols=[]
days=[]

for n in range(num_steps_per_year*years1):
    timeyr1=float(n)/float(num_steps_per_year)
    Sol=1-(math.exp(-timeyr1/timescale1))
    if(n<minilimit): Sol=Sol*minidim
    if(no_dimming==1):  Sol=1
    Sols.append(Sol)
    model1.subprocess['insolation'].S0=S1*Sol
    #if(n%5==0): print(".")
    Ts=model1.Ts
    tmean1=climlab.global_mean(model1.Ts)
    tmeans1.append(tmean1)
    #Tss.append(Ts)
    Tss[:,n] = np.squeeze(Ts)
    days.append(n)
    model1.step_forward()

Sols=np.array(Sols)
days=np.array(days)
years=np.copy(days)/90

Tss=np.array(Tss)[:,:]

#Tss=np.array(Tss)[0,365] ## first year

sel2=snapshot_year*90+1
sel1=sel2-90

Tss2=Tss[0:90,sel1:sel2]

print(np.shape(Tss))

#quit(-1)

Tssmeanprofile=np.sum(Tss, axis=1)

print(np.shape(Tssmeanprofile))

#print(np.shape(Tss))

#quit(-1)

lat1=model1.lat
albedo1=model1.albedo
Ts=model1.Ts

meanTs=climlab.global_mean(model1.Ts)
icelat=np.max(model1.icelat)

print(" Ts_avg ", round(np.mean(meanTs),2))
print(" Icelat ", round(icelat,2))

asr1=model1.ASR
olr1=model1.OLR
albedo1=model1.albedo

print (np.mean((asr1*albedo1)/olr1))

print(np.mean(asr1), np.mean(olr1), np.mean(albedo1))

print (np.shape(Tss))

plt.axhline(y=0, linestyle="--", lw=2)
plt.axhline(y=Tequil_mean, color="green", linestyle=":", lw=2)

plt.plot(years, tmeans1, lw=3, color="blue")
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
#plt.tick_params(axis='both', which='minor', labelsize=15)
plt.title("Impact winter", fontsize=18)
plt.xlabel('Years', fontsize=15)
plt.ylabel('Surface temperature (ºC)', fontsize=15)

plt.show()

#plt.plot(years, tmeans1)
#plt.plot(days,tmeans1)
#plt.plot(years,Sols)

#plot_tann(Ts, lat1)
plotmonths(Tss2, lat1)

plt.show()

latitude=85

plotmonthcurve(Tss2, lat1, latitude)

plt.show()