File:Permian triassic boundary earth mean temperature as function of log10 co2 1.png
Original file (767 × 651 pixels, file size: 62 KB, MIME type: image/png)
Captions
Summary
[edit]DescriptionPermian triassic boundary earth mean temperature as function of log10 co2 1.png |
English: Permian-Triassic boundary Earth - simulated mean temperature as function of log10 CO2 |
Date | |
Source | Own work |
Author | Merikanto |
This image is based on Scotese Paleodems, Exoplasim and various carbon dioxide estimations of Permian-Triassic boundary.
Published August 1, 2018 | Version v2
PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic
Scotese, Christopher R Wright, Nicky M
https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/
https://zenodo.org/records/5460860
Visualization code. Data taken from tas, that is seen on Panoply view of different simulations.
- 3
-
- simulation result tmean, log co2 plot
- 9.12.2023 v 0000.00
-
import numpy as np
import math
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
co2s=np.array([450,1000,3000,5000, 10000]).flatten()
tass=np.array([16.0, 19.3,24.3,26.7,30.7]).flatten()
logco2s=np.log10(co2s)
- model = LinearRegression()
model1 = LinearRegression().fit(logco2s.reshape((-1, 1)), tass.reshape((-1, 1)))
print(f"intercept: {model1.intercept_}")
print(f"slope: {model1.coef_}")
- r_sq = model.score(x, y)
- print(f"coefficient of determination: {r_sq}")
logco2s_new=np.arange(2,5).flatten()
- tass_new = np.arange(0,6).reshape((-1, 1))
- tass_new = np.arange(0,6).reshape((-1, 1))
- tass_new = model.predict(logco2s_new.reshape((-1, 1)))
tass_new = model1.predict(logco2s_new.reshape((-1, 1)) )
plt.scatter(logco2s, tass, s=150, color="red", label="Simulations of P-T")
plt.plot(logco2s_new, tass_new, linestyle="--", color="black" , lw=3, label="Estimated T/CO2 function of P-T")
plt.scatter(math.log10(278), 13.9, s=250, color="green", label= "pre-industrial Earth")
plt.suptitle("Mean temperature of Earth as function of log10 (CO2 ppmv) ", fontsize=15)
plt.title(" Permian-Triassic boundary 251 Ma, Celsius", fontsize=14)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.ylabel("Tmean of Earth (degC)", fontsize=15)
plt.xlabel("log10 CO2 (ppmv) ", fontsize=15)
plt.grid()
plt.legend(fontsize=15)
plt.show()
Exoplasim
https://github.com/alphaparrot/ExoPlaSim
Temperature visualization code is here
https://commons.wikimedia.org/wiki/File:Permian_triassic_boundary_251ma_co2_2500_tas_1.png
Exoplasim code
-
- Exoplasim planet running code, python3, ubuntu
- attempt to create exoplasim restart code
- you can continue running
- based on previous run.
-
- 08.12.2023 0000.0006b
-
- convert to T21, input netcdf
- load one lon, lat, z grid
- or Tarasov glac1d grid
-
- MPI NOTE: if you use more than
-
- one processor, you cannot in most cases run MPI in root
- you can use even number of process in mpi: 2, 4, 6 ..
-
- in ubuntu you must install
-
- pip3 install exoplasim[netCDF4]
- not
- "sudo pip3 install exoplasim[netCDF4]"
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d
import netCDF4
import exoplasim as exo
NLAT=0
NLON=0
def writeSRA(name,kcode,field,NLAT,NLON):
label=name+'_surf_%04d.sra'%kcode
header=[kcode,0,20170927,0,NLON,NLAT,0,0]
fmap = field.reshape((int(NLAT*NLON/8),8))
sheader =
for h in header:
sheader+=" %11d"%h
lines=[]
i=0
while i<NLAT*NLON/8:
l=
for n in fmap[i,:]:
l+=' %9.3f'%n
lines.append(l)
i+=1
text=sheader+'\n'+'\n'.join(lines)+'\n'
f=open(label,'w')
f.write(text)
f.close()
print (label)
def writeSRA2(label,kcode,field,NLAT,NLON):
#label=name+'_surf_%04d.sra'%kcode
header=[kcode,0,20170927,0,NLON,NLAT,0,0]
fmap = field.reshape((int(NLAT*NLON/8),8))
sheader =
for h in header:
sheader+=" %11d"%h
lines=[]
i=0
while i<NLAT*NLON/8:
l=
for n in fmap[i,:]:
l+=' %9.3f'%n
lines.append(l)
i+=1
text=sheader+'\n'+'\n'.join(lines)+'\n'
f=open(label,'w')
f.write(text)
f.close()
print (label)
def savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1):
nlat1=len(xoutlats1)
nlon1=len(xoutlons1)
#indata_set1=indata1
print(outfilename1)
ncout1 = netCDF4.Dataset(outfilename1, 'w', format='NETCDF4')
outlat1 = ncout1.createDimension('lat', nlat1)
outlon1 = ncout1.createDimension('lon', nlon1)
outlats1 = ncout1.createVariable('lat', 'f4', ('lat',))
outlons1 = ncout1.createVariable('lon', 'f4', ('lon',))
outvalue1 = ncout1.createVariable(outvarname1, 'f4', ('lat', 'lon',))
outvalue1.units = 'Unknown'
outlats1[:] = xoutlats1
outlons1[:] = xoutlons1
outvalue1[:, :] =xoutvalue1[:]
ncout1.close()
return 0
def loadnetcdf_single_tomem(infilename1, invarname1):
global cache_lons1
global cache_lats1
print(infilename1)
inc1 = netCDF4.Dataset(infilename1)
inlatname1="lat"
inlonname1="lon"
inlats1=inc1[inlatname1][:]
inlons1=inc1[inlonname1][:]
cache_lons1=inlons1
cache_lats1=inlats1
indata1_set1 = inc1[invarname1][:]
dim1=indata1_set1.shape
nlat1=dim1[0]
nlon1=dim1[1]
inc1.close()
return (indata1_set1)
def create_sras(topo):
global NLAT
global NLON
topo2=np.copy(topo)
masko=np.copy(topo)
topo2[topo2 < 1] = 0
masko[masko < 1] = 0
masko[masko > 0] = 1
grid=np.flipud(masko)
name="Example"
writeSRA(name,129,topo,NLAT,NLON)
writeSRA(name,172,grid,NLAT,NLON)
writeSRA2("topo.sra",129,topo2,NLAT,NLON)
writeSRA2("landmask.sra",172,grid,NLAT,NLON)
return(0)
def convert_to_t21(infilename1, outfilename1):
global NLAT
global NLON
indimx=361
indimy=181
#indimx=360
#indimy=360
## t21 64x32
shapex=64
shapey=32
NLAT=shapex
NLON=shapey
nc = netCDF4.Dataset(infilename1)
inlats=nc['lat'][:]
inlons=nc['lon'][:]
#print(inlats)
#print(inlons)
latlen=len(inlats)
lonlen=len(inlons)
#print(lonlen, latlen)
indimx=lonlen
indimy=latlen
dem=nc['z']
#dem=np.flipud(dem000)
dem2=np.copy(dem)
#dem2[dem2 < 0] = 0
#plt.imshow(dem,cmap='gist_earth')
#plt.imshow(dem2,cmap='gist_earth')
#plt.show()
#quit(0)
lts=[85.7606, 80.2688, 74.7445, 69.2130, 63.6786, 58.1430, 52.6065, 47.0696,
41.5325,35.9951, 30.4576, 24.9199, 19.3822, 13.8445, 8.3067, 2.7689,
-2.7689, -8.3067, -13.8445, -19.3822, -24.9199, -30.4576, -35.9951, -41.5325,
-47.0696, -52.6065, -58.1430, -63.6786, -69.2130, -74.7445, -80.2688, -85.7606]
##
lns=[0, 5.6250, 11.2500, 16.8750, 22.5000, 28.1250, 33.7500 ,39.3750,
45.0000, 50.6250, 56.2500, 61.8750, 67.5000, 73.1250, 78.7500, 84.3750,
90.0000, 95.6250, 101.2500, 106.8750, 112.5000, 118.1250, 123.7500, 129.3750,
135.0000, 140.6250, 146.2500, 151.8750, 157.5000, 163.1250, 168.7500, 174.3750,
180.0000, 185.6250, 191.2500, 196.8750, 202.5000, 208.1250, 213.7500, 219.3750,
225.0000, 230.6250, 236.2500, 241.8750, 247.5000, 253.1250, 258.7500, 264.3750,
270.0000, 275.6250, 281.2500, 286.8750, 292.5000, 298.1250, 303.7500, 309.3750,
315.0000, 320.6250, 326.2500, 331.8750, 337.5000, 343.1250, 348.7500, 354.3750]
ly2=len(lts)
lx2=len(lns)
shapex=lx2
shapey=ly2
#print("sheip")
#print(shapex, shapey)
lons, lats = np.meshgrid(lns,lts)
#print (lts)
#print (lns)
new_W, new_H = (shapey,shapex)
xrange = lambda x: np.linspace(0, 360, x)
f2 = interp2d(xrange(indimx), xrange(indimy), dem2, kind="linear")
#f2 = interp2d(range(indimx), range(indimy), dem2, kind="cubic")
demo = f2(xrange(shapex), xrange(shapey))
#plt.imshow(demo)
#plt.show()
#quit(0)
f3 = interp2d(xrange(indimx), xrange(indimy), dem2, kind="linear")
#masko = f3(xrange(shapex), xrange(shapey))
#topo=np.flipud(demo)
topo=np.copy(demo)
#grid=np.fliplr(masko)
#def savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1):
savenetcdf_single_frommem(outfilename1, "z", topo,lts,lns)
return(topo,lons,lats)
def load_glac1d_dem(indatafile, outdatafile, a_yr):
# load dem from Tarsaov GLAC1d anno domini 2021
global NLAT
global NLON
yr=a_yr
lok=int(abs(yr/100-260))
# tarasov ice 26k
nc = netCDF4.Dataset(indatafile1)
#print(nc)
eisbase=nc['ICEM']
inlats=nc['YLATGLOBP5'][:]
inlons=nc['XLONGLOB1'][:]
dem=nc['HDCB'][lok]
#dem=np.flipud(dem000)
#print (dem)
#print (np.shape(dem))
#plt.imshow(dem,cmap='gist_earth')
savenetcdf_single_frommem(outdatafile, "z",dem,inlats,inlons)
return(0)
- maybe nok
def convert_to_t42(infilename1, outfilename1):
## ONLY attempi! to create T42!
global NLAT
global NLON
indimx=361
indimy=181
## t42 64x32
#shapex=64
#shapey=32
shapex=128
shapey=64
#shapey=63
NLAT=shapex
NLON=shapey
nc = netCDF4.Dataset(infilename1)
inlats=nc['lat'][:]
inlons=nc['lon'][:]
latlen=len(inlats)
lonlen=len(inlons)
indimx=lonlen
indimy=latlen
dem=nc['z']
#dem=np.flipud(dem000)
dem2=np.copy(dem)
## test t21
tdx=360.0/shapex
#tdy=180.0/shapey
tdy=(90.0-85.706)/2
minix=0.0
maksix=360-tdx
maksiy=90-tdy
miniy=-90+tdy
#print(90-tdy)
#
#print(miniy)
#print(maksiy)
#quit(-1)
#lns=np.linspace(minix, maksix, num=shapex)
#lts=np.linspace(maksiy, miniy, num=shapey)
## jn WARNING 90!
lts=[87.8638, 85.0965 ,82.3129, 79.5256, 76.7369 ,73.9475 ,71.1578, 68.3678, #ok
65.5776, 62.7874, 59.9970 ,57.2066, 54.4162, 51.6257, 48.8352, 46.0447,
43.2542, 40.4636, 37.6731 ,34.8825, 32.0919, 29.3014, 26.5108, 23.7202,
20.9296, 18.1390, 15.3484 ,12.5578, 9.7671, 6.9765, 4.1859, 1.3953,
-1.3953, -4.1859, -6.9765, -9.7671, -12.5578, -15.3484, -18.1390, -20.9296,
-23.7202,-26.5108, -29.3014 ,-32.0919, -34.8825, -37.6731, -40.4636,-43.2542,
-46.0447,-48.8352, -51.6257, -54.4162, -57.2066, -59.9970, -62.7874, -65.5776,
-68.3678,-71.1578 ,-73.9475, -76.7369 ,-79.5256, -82.3129, -85.0965, -87.8638]
lns=[0.0000 ,2.8125, 5.6250, 8.4375, 11.2500, 14.0625 ,16.8750 ,19.6875,
22.5000,25.3125, 28.1250, 30.9375 ,33.7500,36.5625 ,39.3750, 42.1875,
45.0000,47.8125, 50.6250, 53.4375, 56.2500, 59.0625 ,61.8750, 64.6875,
67.5000, 70.3125, 73.1250, 75.9375, 78.7500, 81.5625, 84.3750, 87.1875,
90.0000, 92.8125, 95.6250 ,98.4375 ,101.2500, 104.0625, 106.8750, 109.6875,
112.5000, 115.3125, 118.1250, 120.9375,123.7500 ,126.5625 ,129.3750, 132.1875,
135.0000, 137.8125, 140.6250 ,143.4375, 146.2500 ,149.0625, 151.8750 ,154.6875,
157.5000, 160.3125, 163.1250, 165.9375, 168.7500, 171.5625 ,174.3750, 177.1875,
180.0000, 182.8125, 185.6250 ,188.4375, 191.2500, 194.0625, 196.8750, 199.6875,
202.5000, 205.3125, 208.1250, 210.9375, 213.7500 ,216.5625, 219.3750 ,222.1875,
225.0000, 227.8125, 230.6250 ,233.4375, 236.2500, 239.0625, 241.8750, 244.6875,
247.5000, 250.3125, 253.1250, 255.9375, 258.7500, 261.5625, 264.3750, 267.1875,
270.0000, 272.8125, 275.6250, 278.4375, 281.2500 ,284.0625 ,286.8750, 289.6875,
292.5000, 295.3125, 298.1250, 300.9375, 303.7500 ,306.5625, 309.3750, 312.1875,
315.0000, 317.8125, 320.6250, 323.4375, 326.2500, 329.0625 ,331.8750, 334.6875,
337.5000, 340.3125, 343.1250, 345.9375, 348.7500, 351.5625 ,354.3750 ,357.1875]
#lns=
#print (lts)
#print (lns)
#print (len(lns),len(lts))
#quit(-1)
ly2=len(lts)
lx2=len(lns)
shapex=lx2
shapey=ly2
#print("sheip")
#print(shapex, shapey)
lons, lats = np.meshgrid(lns,lts)
new_W, new_H = (shapey,shapex)
xrange = lambda x: np.linspace(0, 360, x)
f2 = interp2d(xrange(indimx), xrange(indimy), dem2, kind="linear")
demo = f2(xrange(shapex), xrange(shapey))
f3 = interp2d(xrange(indimx), xrange(indimy), dem2, kind="linear")
topo=demo
savenetcdf_single_frommem(outfilename1, "z", topo,lts,lns)
return(topo,lons,lats)
- exoplasim ,,,
def exo_runner_restarting(firstrun,a_input_dem1, a_gridtype, a_layers, a_years,a_timestep,a_snapshots,a_ncpus,a_eccentricity,a_obliquity,a_lonvernaleq,a_pCO2):
output_format=".nc"
a_pO2=1-a_pCO2-0.79
a_pN2=(1-0.21-a_pCO2)
print("Process input grid, to type ",a_gridtype)
if(a_gridtype=="T21"):
print("T21")
topo, lons, lats=convert_to_t21(a_input_dem1,"demT21.nc")
if(a_gridtype=="T42"):
print("T42")
topo, lons, lats=convert_to_t42(a_input_dem1, "demT42.nc")
create_sras(topo)
print("Creating exoplasim object ")
testplanet= exo.Earthlike(workdir="planet_run",modelname="PLANET",ncpus=a_ncpus,resolution=a_gridtype,layers=a_layers, outputtype=output_format, crashtolerant=True)
glaciers1= {
"toggle": True,
"mindepth":2,
"initialh":-1
}
fluxi1=1338
testplanet.configure(
startemp=5772.0,
flux=fluxi1,# Stellar parameters
eccentricity=a_eccentricity,
obliquity=a_obliquity,
lonvernaleq=a_lonvernaleq,
fixedorbit=True, # Orbital parameters
rotationperiod=1, # Rotation
topomap="topo.sra",
landmap="landmask.sra",
radius=1.0,
gravity=9.80665,
#stormclim=False,
vegetation=2, #toggles vegetation module; 1 for static vegetation, 2 to allow growth
vegaccel=1,
seaice=True,
maxsnow=-1,
glaciers=glaciers1,
pN2=a_pN2,
pCO2=a_pCO2,
pO2=a_pO2,
ozone=True, # Atmosphere
timestep=a_timestep,
snapshots=0, ## jos a_snapshots, vie muistia!
wetsoil=True,
physicsfilter="gp|exp|sp",
restartfile="ressus"
) # Model dynamics
testplanet.exportcfg()
runc1=1
n=0
if(firstrun==1):
print("Creating first restart.")
print("Running ExoPlasim ... ")
testplanet.run(years=1,crashifbroken=True)
lon = testplanet.inspect("lon")
lat = testplanet.inspect("lat")
ts =testplanet.inspect("tsa",tavg=True)
tsavg=np.mean(ts)-273.15
print("Year: ",n," tsa: ",tsavg)
savename = 'ressu'
testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True)
testplanet.save(savename)
looplen=a_years1
peen=0
runc1=1
for n in range(0,looplen):
print("Loop year ",n)
testplanet.modify(flux=fluxi1) #number of output times (months) in the output files
testplanet.exportcfg()
runc1=1
testplanet.run(years=1,crashifbroken=True)
lon = testplanet.inspect("lon")
lat = testplanet.inspect("lat")
ts =testplanet.inspect("tsa",tavg=True)
tsavg=np.mean(ts)-273.15
print("Year: ",n," tsa: ",tsavg)
savename = 'ressu'+str(runc1)
testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True)
testplanet.save(savename)
print("Return.")
return(0)
-
-
print(" Exoplasim simulation restart code ---")
- jn warning maybe nok
- input_dem='./indata/indem.nc'
- input_dem='./indata/Map22_PALEOMAP_1deg_Mid-Cretaceous_95Ma.nc'
- input_dem='./indata/Map14_PALEOMAP_1deg_Paleocene_Eocene_Boundary_55Ma.nc'
- input_dem='/indata/Map13_PALEOMAP_1deg_Early_Eocene_50Ma.nc'
- input_dem='./indata/Map12_PALEOMAP_1deg_early_Middle_Eocene_45Ma.nc'
- input_dem='./indata/Map18_PALEOMAP_1deg_Late_Cretaceous_75Ma.nc' ## OK
- input_dem='./indata/Map20_PALEOMAP_1deg_Late_Cretaceous_85Ma.nc' ## nok
- input_dem='./indata/Map24_PALEOMAP_1deg_Early Cretaceous_105Ma.nc' ## nok
- input_dem='./indata/Map17_PALEOMAP_1deg_Late_Cretaceous_70Ma.nc' ##nok
- input_dem='./indata/Map19_PALEOMAP_1deg_Late_Cretaceous_80Ma.nc'
- input_dem="./indata/Map16_PALEOMAP_1deg_KT_Boundary_65Ma.nc"
- input_dem="./indata/Map43_PALEOMAP_1deg_Late_Triassic_200Ma.nc"
- input_dem='./indata/Map19_PALEOMAP_1deg_Late_Cretaceous_80Ma.nc' ## OK
- input_dem='./indata/Map21_PALEOMAP_1deg_Mid-Cretaceous_90Ma.nc' #90ma
input_dem='./maps1/Map49_PALEOMAP_1deg_Permo-Triassic Boundary_250Ma.nc' # PT raja co2 1600. jopa 3000-4000
- input_dem='./indata/Map57_PALEOMAP_1deg_Late_Pennsylvanian_300Ma.nc' ## Late Pennsylcanian ice, co2 200? 250?
- input_dem="./indata/Map56_PALEOMAP_1deg_Early_Permian_295Ma.nc"
- indatafile1='./indata/TOPicemsk.GLACD26kN9894GE90227A6005GGrBgic.nc'
- input_dem="origodem.nc"
- a_yr=14500
- load_glac1d_dem(indatafile1, input_dem, 14500)
- input one de scotese palaeomap dem!
- def convert_to_t42(infilename1, outfilename1):
- topo, lons, lats=convert_to_t21(input_dem, "demT21.nc")
- topo, lons, lats=convert_to_t42(input_dem, "demT42.nc")
- plt.imshow(topo,cmap='gist_earth')
- plt.show()
- input_dem="./sand.nc" ##dem of desert planet
a_modelname1="planet"
a_workdir1="planet_run"
a_runsteps1=200
a_years1=a_runsteps1
a_timestep1=30
a_snapshots1=0
a_ncpus1=4
a_layers1=8
a_outputtype1=".nc"
- a_resolution1="T42"
a_resolution1="T21"
a_precision1=4
a_crashtolerant1=True
a_landmap1="landmask.sra"
a_topomap1="topo.sra"
- nowadays ca 0 BP
- a_eccentricity1=0.01671022
- a_obliquity1=23.44
- a_lonvernaleq1=102.7
- a_pCO21=360e-6
- 10000 yrs ago
- a_eccentricity1=0.0194246086670259
- a_obliquity1=24.230720588
- a_lonvernaleq1=295.26651297
- a_pCO21=265e-6
- 14500 yrs ago
- a_eccentricity1=0.019595
- a_obliquity1=23.6801
- a_lonvernaleq1=221.5
- (229.64+213.3)/2
- a_pCO21=210e-6
- 25000 yrs ago
- a_eccentricity1=0.0178681374211005
- a_obliquity1= 22.408850897
- a_lonvernaleq1=49.92
- a_pCO21=180e-6
- permo
a_eccentricity1=0.0167022
a_obliquity1=23.441
a_lonvernaleq1=282.7
a_pCO21=450.0e-6
- a_pCO21=700.0e-6
- early permian 295 ma
- late pennsylvanian 300 ma
- a_eccentricity1=0.01671022
- a_obliquity1=23.441
- a_lonvernaleq1=282.7
- a_pCO2=250.0e-6 ## ca 200 - 250 ppmvol
- a_pCO21=180.0e-6
- a_pCO21=100.0e-6
- permo-triassic boundary ca 250 ma
- a_eccentricity1=0.01671022
- a_obliquity1=23.441
- a_lonvernaleq1=282.7
- a_pCO21=1600.0e-6 ## cal1600 ppmvol 3000 ? 2000-4000
print("Exoplasim ...")
- if you run simu first time, you must set
- firstrun=1
- firstrun=1
firstrun=1
a_years1=500
exo_runner_restarting(firstrun, input_dem, a_resolution1, a_layers1, a_years1,a_timestep1,a_snapshots1,a_ncpus1,a_eccentricity1,a_obliquity1,a_lonvernaleq1,a_pCO21)
print(".")
Licensing
[edit]- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 18:35, 9 December 2023 | 767 × 651 (62 KB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Software used | |
---|---|
Horizontal resolution | 39.37 dpc |
Vertical resolution | 39.37 dpc |