File:Temperature of fast rotaing desert planet by latitude simple ebm and diffusion 1 1 1 1.png

https://brian-rose.github.io/ClimateLaboratoryBook/courseware/numerical-diffusion.html

1. -*- coding: utf-8 -*

1. 2. fast rotating desert planet temperature vs latitude
   3. simple basic energy balance temperature
2. 1. insolation effect, albedo, diffusion
3. 2. python 3 source code
4. 1. 03.11.2023 v 0000.0000

import numpy as np import scipy

import math import matplotlib.pyplot as plt import PIL

from scipy.ndimage import convolve from scipy.integrate import odeint

def diffusive_flux(u, deltax, K=1): # Take the finite difference F = np.diff(u)/deltax # add a zero as the first element (no flux on boundary) F = np.insert(F, 0, 0.) # add another zero as the last element (no flux on boundary) F = np.append(F, 0.) # flux is DOWN gradient, proportional to D return -K*F

def diffusion(u, deltax, K=1): # compute flux F = diffusive_flux(u, deltax, K) # take convergence of flux return -np.diff(F) / deltax

def step_forward(u, deltax, deltat, K=1): dudt = diffusion(u, deltax, K) return u + deltat * dudt

def calc_slowrot_planet_water_line(temp_K): theta=math.acos( math.pow((temp_K/298),4) )*(180/3.141828) return(theta)

1. 1. main code

numlats=180

Kelvin=273.15

albedo=0.3

1. albedo=0.3 ## earth

sb=5.670374419e-8

1. S1coeff=0.65

S1coeff=1.0

S1=1361.5*S1coeff S14=S1/4

1. S1=1365.2*1.0

1. print(1/math.sqrt(S1coeff))

lats=np.linspace(-90,90,numlats) latsrad=lats*math.pi/180 latsrad0=np.copy(latsrad)

cosun=np.copy(latsrad)

cosun=np.where(cosun<(-math.pi/2),(-math.pi/2), cosun ) cosun=np.where(cosun>(math.pi/2),(math.pi/2), cosun )

zenitarads=abs(np.copy(cosun))

1. zenitarads=abs(math.pi/2-abs(zenitarads))

cosun=np.cos(cosun)

1. plt.plot(lats,zenitarads)

1. plt.show()

1. quit(1)

no_absorb_scatter=0.925

insols=cosun*S1*no_absorb_scatter

1. plt.plot(lats,insols)
2. plt.plot(lats,insols2)

1. plt.show()

1. quit(1)

energys0=np.copy(insols) energys1=np.copy(energys0)

teqs_1=np.power(((energys1*(1-albedo)) / (4*sb)), 1/4)

teq_planet_1=math.pow(((S1*(1-albedo)) / (4*sb)), 1/4)

1. 1. https://en.wikipedia.org/wiki/Idealized_greenhouse_model

1. 1. greenhouse tsurf

tss_1=1.189*teqs_1

teqs=teqs_1

1. iceline_slow_1=calc_slowrot_planet_water_line(teq)
2. tem1=np.mean(teqs)

temps0=np.copy(tss_1) temps1=np.copy(tss_1)

tems1=np.mean(temps0)+temps0*0

degrad1=math.pi/180

1. for n in range(0, 180):

tempscache0=np.copy(temps1) tempscache1=np.copy(temps1)

u=np.copy(tempscache0)

1. 1. u, dx, dt , k

dx=1 dt=1 K=0.3 ## for dx=1 deg, dt=1, 1000 steps

for n in range(0,1000): u=step_forward(u, dx,dt,K) #print(np.min(u),np.max(u))

tempscache1=np.copy(u)

temps0=np.copy(temps1) temps1=np.copy(tempscache1)

1. plt.plot(tempscache0)

1. plt.plot(tempscache1)
2. plt.show()

tdmax1=round((np.max(temps1-Kelvin)),1) tdmin1=round((np.min(temps1-Kelvin)),1)

temps_cache=np.copy(temps1) latsrad=np.copy(latsrad0) coslats=np.cos(latsrad) coslats2=np.power(coslats,2) temps_cache=temps_cache-Kelvin

temps_cache=(temps_cache)*np.abs(coslats2) tdmean1=round((np.mean(temps_cache)),1) tdelta1=round((abs(tdmax1-tdmin1)),1)

1. print(coslats)

print(round((tdmean1),1),round((tdmax1),1), round((tdmin1),1), round((tdelta1),1))

1. quit(-1)

idxts0 = np.where(temps1 >(0+Kelvin) )

1. print(idxts0)

idxt01=np.array(list(idxts0[0]))

1. print (idxt01)

idxta0=idxt01[0] idxtb0=idxt01[-1]

1. print(idxta0)

latta0= round(lats[idxta0],1) lattb0= round(lats[idxtb0],1)

1. print("Icelats ",latta0, lattb0)

iceline_1=latta0

1. print(tem1)

1. quit(-1)

print("Iceline lat deg " ,round(iceline_1,1) )

1. quit(-1)

title1="Temperature of fast rotating desert planet \n " title2="°C if S="+ str(round(S1,2))+ " Wm², a="+str(round(albedo,3))+" , greenhouse "

1. +" "+ str(tdmean1)," ",str(tdmax1)," ", str(tdmin1)," ", str(tdelta1)

plt.xlabel("latitude", size=16) plt.ylabel("Temperature degrees Celsius", size=16) plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.title(title2, size=12, color="#003f00")

plt.suptitle(title1, size=18)

tmaxlabel1="Tmax "+str(round(tdmax1-Kelvin,1))+ " °C"

icelatlabel1="Icelat " +str(round(lattb0,1))+ " °"

plt.xlim(-100,100) plt.hlines(xmin=-180, xmax=180, y=0, color="blue",lw=2, linestyle=":", label="Ice/water")

1. plt.hlines(xmin=-180, xmax=180, y=-10, color="blue", linestyle=":")
   1. co2 ice

plt.hlines(xmin=-180, xmax=180, y=-78, color="lightblue",lw=2, linestyle=":", label="CO² freezes") plt.plot(lats,temps1-273.15, lw=5, color="orange",linestyle="-", label="Surface temperature °C")

plt.scatter(0,(tdmax1),s=300, marker="o", color="red", label=tmaxlabel1 )

1. plt.vlines(ymin=-40, ymax=0, x=latt0, color="blue",lw=3, linestyle="--")

plt.scatter(latta0, 0, s=300,marker="X", color="blue")

plt.scatter(lattb0, 0, s=300,marker="X",color="blue",label=icelatlabel1)

1. plt.text(-90,tdmax1,"Tmax"+str(tdmax1))

1. plt.plot(lats, teqs_1-Kelvin)
2. plt.plot(lats, teqs_2-Kelvin)
3. plt.plot(lats, teqs_3-Kelvin)
4. plt.plot(lats, tdesert1-Kelvin)

plt.grid() plt.legend(fontsize=10) plt.show()