File:Permian triassic boundary 251ma co2 2500 tas 1.png

Scotese, Christopher R  Wright, Nicky M
 

2. 1. temperature plotter
3. python 3
   1. 09.12.2023 0000.0006a

import matplotlib.mlab as mlab import matplotlib.pyplot as plt import netCDF4 from netCDF4 import Dataset from matplotlib.pylab import * import numpy as np from scipy import interpolate from scipy.interpolate import griddata

from sklearn.linear_model import LinearRegression from pygam import LogisticGAM from pygam import LinearGAM from sklearn.ensemble import RandomForestRegressor from sklearn.preprocessing import StandardScaler from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import make_pipeline

1. from skimage.transform import resize

import skimage

1. 1. random forest setting

RF_estimators1=1 RF_features1=1

def gam_multiple_vars( x1, y1, x2): print ("GAM") xa1=np.array(x1) ya1=np.array(y1) xa2=np.array(x2)

#print (xa1)

1. quit(-1)

sha1=np.shape(x1)

dim2=sha1[1]

x=xa1.reshape((-1, dim2)) y=ya1.reshape((-1, 1)) xb=xa2.reshape((-1, dim2))

#sc = StandardScaler() #xa = sc.fit_transform(x) #xba = sc.transform(xb)

#model = LogisticGAM().fit(x, y) model = LinearGAM().fit(x, y)

y2= model.predict(xb)

return(y2)

def random_forest_multiple_vars( x1, y1, x2): print ("RF")

global RF_estimators1 global RF_features1

print(RF_estimators1, RF_features1)

#quit(-1)

xa1=np.array(x1) ya1=np.array(y1) xa2=np.array(x2)

#print (xa1)

1. quit(-1)

sha1=np.shape(x1)

dim2=sha1[1]

x=xa1.reshape((-1, dim2)) y=ya1.reshape((-1, 1)) xb=xa2.reshape((-1, dim2))

#model = LinearRegression().fit(x, y) #degree=3 #polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression()) #model=polyreg.fit(x,y)

sc = StandardScaler() xa = sc.fit_transform(x) xba = sc.transform(xb)

## orig ##model = RandomForestRegressor(n_estimators=10, max_features=2).fit(xa,y) model = RandomForestRegressor(n_estimators=RF_estimators1, max_features=RF_features1).fit(xa,y)

y2= model.predict(xba) return(y2)

def linear_regression_multiple_vars( x1, y1, x2): print ("MLR") xa1=np.array(x1) ya1=np.array(y1) xa2=np.array(x2)

sha1=np.shape(x1)

dim2=sha1[1]

x=xa1.reshape((-1, dim2)) y=ya1.reshape((-1, 1)) xb=xa2.reshape((-1, dim2))

#model = LinearRegression().fit(x, y) degree=3 polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression()) model=polyreg.fit(x,y)

y2= model.predict(xb) return(y2)

def linear_regression_singe_var( x1, y1, x2): #print (x1) #print(y1) #return(0)

#xa1=np.array(x1) #ya1=np.array(y1) #xa2=np.array(x2) xa1=np.array( x1) ya1=np.array(y1) xa2=np.array(x2)

sha1=np.shape(x1)

dim2=sha1[1]

x=xa1.reshape((-1, dim2)) y=ya1.reshape((-1, 1)) xb=xa2.reshape((-1, dim2)) #x=xa1.reshape((-1, 1)) #y=ya1.reshape((-1, 1)) #xb=xa2.reshape((-1, 1))

#print (x) #print (y)

model = LinearRegression().fit(x, y) #model = RandomForestRegressor(n_estimators=10, max_features=2).fit(x,y) #degree=2 #polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression()) #polyreg=make_pipeline(PolynomialFeatures(degree), ) #model=polyreg.fit(x,y)

## warning not xb y2= model.predict(xb) #print(y2) #print("LR") return(y2)

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="latitude" inlonname1="longitude" 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)

2. main

1. kaption1="Temperature of July 67 million years ago, degrees Celsius"

1. kaption2="S=1312 CO2=180 ppm obliq=22 ecc=0.06 mvelp=90"

a_eccentricity1=0.0167022 a_obliquity1=23.441 a_lonvernaleq1=282.7

1. a_pCO21=5000.0e-6

kaption1="Mean temperature of year degrees Celsius, PTME 251 Ma " kaption2="S=1338 CO2=2500 ppm obliq=23.441 ecc=0.0167 mvelp=282.7"

fintemp1="./simut1/pt_250_co2_2500_mvelp_282.nc"

outname1="out.png"

1. outname1="cretaceous_67ma_co2_4x_tjanuary_1.png"

year1=0 month1=0 ## if 0, mean surf temp of year

1. month1=1

vari1="tas"

tsa_offset=0

1. fintemp1="./exoplasim_445ma_co2_125/MOST.00040.nc"

1. findem1="Map22_PALEOMAP_6min_Mid-Cretaceous_95Ma.nc"

findem1="./maps1/Map49_PALEOMAP_6min_Permo-Triassic Boundary_250Ma.nc"

1. findem1="./maps/Map49_PALEOMAP_1deg_Permo-Triassic Boundary_250Ma.nc"

1. 1. downscaleing size

dimy3=180 dimx3=360

extent1=(-180,180,-90,90)

1. extent1=(-180,180,90,-90)

1. print(num1, kk1)

1. print(fintemp1)

kindex1=1

ncinu1 = Dataset(fintemp1, "r+", format = "NETCDF4")

1. print(ncinu1)

1. quit(-1)

ncdem1 = Dataset(findem1, "r+", format = "NETCDF4")

1. print(ncdem1)

1. quit(-1)

1. navalue=9.969209968386869e+36

lon0 = ncinu1.variables['lon'] lat0 = ncinu1.variables['lat'] lon1 = ncdem1.variables['longitude'] lat1 = ncdem1.variables['latitude']

dimx1=len(lon0) dimy1=len(lat0)

dimx2=len(lon1) dimy2=len(lat1)

tsa00 = np.array(ncinu1.variables[vari1])

1. snd = np.array(ncinu1.variables['snd'])

tsa=np.mean(tsa00, axis=0)

1. if(month1==0): tsa=np.mean(tsa, axis=0)
2. else:tsa=tsa[month1-1]

tsa=tsa-273.15

tsa_january=np.array(ncinu1.variables['tsa'][0])-273.15 tsa_july=np.array(ncinu1.variables['tsa'][6])-273.15

tsa_min0=round(np.min(tsa),2) tsa_max0=round(np.max(tsa),2)

1. tsa_mean=np.mean(tsa)

coslat0=np.cos(np.radians(lat0))

tsalat0=np.mean(tsa, axis=1) print (tsalat0)

costemp0=coslat0*coslat0*tsalat0

1. tsa_mean0=round(np.mean(costemp0),2)

tsa_mean0=round(np.mean(tsalat0),2)

print(tsa_mean0)

1. quit(-1)

mask00=ncinu1.variables['lsm'][1] mask1=np.array(mask00)

snowdepths=np.array(ncinu1.variables['snd'][:]) minsnowdepth=np.min(snowdepths, axis=0) maxsnowdepth=np.max(snowdepths, axis=0)

1. plt.imshow(maxsnowdepth)

1. plt.contour(minsnowdepth, levels=[0])

1. plt.imshow(mask00, cmap="terrain")

1. plt.contour(minsnowdepth, levels=[0])

1. plt.imshow(tsa_january, cmap="jet")

1. plt.imshow(tsa_july, cmap="jet")
2. plt.imshow(tsa, cmap="jet")

1. plt.show()

1. quit(-1)

dem00=np.array(ncdem1.variables['z'])

dem01=np.copy(dem00)

dem1=skimage.transform.resize(dem01,(dimy3, dimx3))

1. 1. warning
   2. dem1=np.flipud(dem1)

1. plt.imshow(np.array(ncinu1.variables['tso'][1]), extent=extent1)

1. plt.imshow(tsa_january, cmap="terrain")

1. plt.imshow(dem1, cmap="terrain")

1. plt.show()

1. quit(-1)

1. print(tsa00)

smalons=np.zeros((dimy1, dimx1))*255 smalats=np.zeros((dimy1, dimx1))*255

for iy in range(dimy1): for ix in range(dimx1): sx=ix/dimx1 sy=iy/dimy1 smalons[iy, ix]=sx smalats[iy, ix]=sy

biglons=np.zeros((dimy3, dimx3))*255 biglats=np.zeros((dimy3, dimx3))*255

for iy in range(dimy3): for ix in range(dimx3): sx=ix/dimx3 sy=iy/dimy3 biglons[iy, ix]=sx biglats[iy, ix]=sy

1. plt.imshow(smalons)
2. plt.imshow(smalats)

1. plt.imshow(biglons)
2. plt.imshow(biglats)

1. plt.show()

1. 1. resize input dem

1. print (len(dem00))

kmap2 = plt.cm.get_cmap('RdBu_r') kmap1=plt.get_cmap('YlGnBu')

dmap1=plt.get_cmap('binary')

1. kmap2 = plt.cm.get_cmap('bwr')

kmap2 = plt.cm.get_cmap('coolwarm') ## ok

1. kmap2 = plt.cm.get_cmap('jet')

1. kmap2 = plt.cm.get_cmap('turbo') ##ok
2. kmap2 = plt.cm.get_cmap('gist_rainbow_r')

lon1=np.array(lon0) lat1=np.array(lat0)

1. print(lon1)

1. tsa=np.array(tsa00)
2. mask1=np.array(mask00) ## ok

1. print(np.shape(tsa[month1]))

1. tsa=tsa+tsa_offset

1. plt.imshow(tsa, cmap=kmap2, interpolation="bicubic")
2. plt.imshow(dem00, cmap=kmap2, interpolation="bicubic")

1. plt.show()

1. quit(-1)

dem2=np.copy(dem1)

dem2[dem2>-1]=1 dem2[dem2<0]=0

1. sla=[]

1. sla.append(tsa)

1. apoints1=
2. bpoints1=list(zip(*sla))

dem4=np.copy(dem1) dem4[dem4<0]=0

masco1=np.copy(dem1) masco1[masco1<0]=0 masco1[masco1>0]=255

1. plt.imshow(masco1)
2. plt.imshow(mask1)

1. plt.imshow(tsa)

1. plt.imshow(dem4)
2. plt.show()

1. quit(-1)

smalldem = skimage.transform.resize(dem4,(dimy1, dimx1)) bigsmalldem=skimage.transform.resize(smalldem,(dimy3, dimx3), anti_aliasing=True)

demdifference=bigsmalldem-dem4

tempdifference=demdifference*10/1000 ## simple assumption lapse rate 10 or 6.5 / 1000 m

1. 1. add some random noise

mu, sigma = 0, 0.3 # mean and standard deviation

random0 = np.random.normal(mu, sigma, (int(dimy3/10), int(dimx3/10))) random1=skimage.transform.resize(random0, (dimy3, dimx3), anti_aliasing=True ) mu1, sigma1 = 0, 0.1 random0 = np.random.normal(mu1, sigma1, (int(dimy3/3), int(dimx3/3))) random2=skimage.transform.resize(random0, (dimy3, dimx3), anti_aliasing=True )

1. 1. tempdifference=tempdifference+random1+random2

bigtsa=skimage.transform.resize(tsa, (dimy3, dimx3), anti_aliasing=True )

dstsa=bigtsa+tempdifference

machinelearn=0

if(machinelearn==1): sla=[] sla.append(smalldem) sla.append(smalats) sla.append(smalons) slb=[] slb.append(dem4) slb.append(biglats) slb.append(biglons) apoints1=list(zip(*sla)) ## check this! bpoints1=tsa cpoints1=list(zip(*slb)) dows=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1) #dows=gam_multiple_vars(apoints1, bpoints1, cpoints1) #dows=random_forest_multiple_vars(apoints1, bpoints1, cpoints1) #print (shape(dows)) dows2=np.reshape(dows, (dimy3, dimx3)) #dows3=dows2-np.min(dows2) #lapses=(dows3/dem4)*1000 #print (np.mean(lapses)) #plt.imshow(dows2) #plt.show() #quit(-1)

1. plt.imshow(mask1)
2. plt.imshow(bigtsa)
3. plt.imshow(bigsmalldem)
4. plt.imshow(bigtsa, cmap=kmap2)

1. plt.imshow(smalldem)
2. plt.imshow(bigsmalldem)

1. plt.imshow(dstsa, cmap=kmap2)

1. plt.imshow(tempdifference, cmap=kmap2)

1. plt.show()

1. print(mask1)

1. quit(-1)

1. avgtemp1=np.mean(tsa)

1. print("Tavg: ",avgtemp1)

1. quit(-1)

newx=64 newy=32

X = np.arange(-180, 180, 360/newx) Y = np.arange(-90, 90, 180/newy)

oname1="atm_temp.nc" ovar1="Band1"

savenetcdf_single_frommem(oname1, ovar1, tsa,Y, X)

X2 = np.arange(-180, 180, 360/dimx3) Y2 = np.arange(-90, 90, 180/dimy3)

oname1="tempds.nc" ovar1="Band1"

savenetcdf_single_frommem(oname1, ovar1, np.flipud(dstsa),Y2, X2)

1. contourrange1=[-50,-30,-20,-10,0,10,20,30,40,50]

contourrange1=[-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50]

contourrange2=[-50,0,50]

contourrange3=np.arange(-10,40,1)

1. range3=np.arange(-50,50,0.1)

range3=np.arange(-10,40,0.1)

koloro1=['#0000FF', '#A0A0A0', '#FF0000']

imsiz1=(16,8)

fig,ax = plt.subplots(figsize=imsiz1)

size1=24

1. kaption3=" Tmin: "+str(tsa_min0)+" Tmax: "+str(tsa_max0)+" Tmean: "+str(tsa_mean0)+ " °C"

kaption3=" Tmin: "+str(tsa_min0)+" Tmax: "+str(tsa_max0)

plt.suptitle(kaption1, fontsize=size1) ax.set_title(kaption2, fontsize=size1-3) size2=21

ax.annotate(kaption3, xy = (0.02, 0.9), xycoords='axes fraction', fontsize=22, alpha=0.5)

ax.set_xlabel("Longitude",fontsize=size2 ) ax.set_ylabel("Latitude", fontsize=size2 ) ax.tick_params(axis='x', labelsize= size2) ax.tick_params(axis='y', labelsize= size2)

tsa=np.flipud(tsa)

dem1=np.flipud(dem1)

mask1=np.flipud(mask1)

dstsa=np.flipud(dstsa) masco1=np.flipud(masco1)

1. plt.margins(0.0015, tight=True)

1. var1=ncinu1.variables['tso'][6]

1. plt.imshow(dstsa, extent=extent1,interpolation="bicubic", cmap="coolwarm")

1. plt.show()
2. quit()

1. selectvar=tsa

selectvar=dstsa

1. selectvar=tsa

ax.imshow( selectvar, origin="lower", extent=extent1, cmap=kmap2, vmin=-50, vmax=60 ,interpolation="bicubic")

cs = ax.contour(selectvar, origin="lower",extent=extent1, inline=True, colors='#7f0000', alpha=0.5, levels=contourrange1)

ax.clabel(cs, fontsize=20,fmt = '%3.0f', colors='#7f0000')

csdem2 = ax.contour(masco1, origin="lower",extent=extent1, linestyles='-',alpha=0.75, linewidths=[4], inline=True, colors=["#003f00"], levels=[-15000,0,15000])

fig.savefig(outname1)

print(tsa_min0, tsa_max0, tsa_mean0)

print("plot")

plt.show()

print(".")

2. 1. Exoplasim planet running code, python3, ubuntu
4. attempt to create exoplasim restart code
   1. you can continue running
   2. based on previous run.
5. 1. 16.06.2022 0000.0006
8. 1. convert to T21, input netcdf
   2. load one lon, lat, z grid
   3. or Tarasov glac1d grid
10. 1. MPI NOTE: if you use more than
11. 1. one processor, you cannot in most cases run MPI in root
    2. you can use even number of process in mpi: 2, 4, 6 ..
12. 1. in ubuntu you must install
13. 1. pip3 install exoplasim[netCDF4]
    2. not
    3. "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)

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

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

1. 1. jn warning maybe nok

1. input_dem='./indata/indem.nc'
2. input_dem='./indata/Map22_PALEOMAP_1deg_Mid-Cretaceous_95Ma.nc'
3. input_dem='./indata/Map14_PALEOMAP_1deg_Paleocene_Eocene_Boundary_55Ma.nc'
4. input_dem='/indata/Map13_PALEOMAP_1deg_Early_Eocene_50Ma.nc'
5. input_dem='./indata/Map12_PALEOMAP_1deg_early_Middle_Eocene_45Ma.nc'
6. input_dem='./indata/Map18_PALEOMAP_1deg_Late_Cretaceous_75Ma.nc' ## OK
7. input_dem='./indata/Map20_PALEOMAP_1deg_Late_Cretaceous_85Ma.nc' ## nok
8. input_dem='./indata/Map24_PALEOMAP_1deg_Early Cretaceous_105Ma.nc' ## nok
9. input_dem='./indata/Map17_PALEOMAP_1deg_Late_Cretaceous_70Ma.nc' ##nok
   1. input_dem='./indata/Map19_PALEOMAP_1deg_Late_Cretaceous_80Ma.nc'

1. input_dem="./indata/Map16_PALEOMAP_1deg_KT_Boundary_65Ma.nc"

1. input_dem="./indata/Map43_PALEOMAP_1deg_Late_Triassic_200Ma.nc"

1. input_dem='./indata/Map19_PALEOMAP_1deg_Late_Cretaceous_80Ma.nc' ## OK

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

1. input_dem='./indata/Map57_PALEOMAP_1deg_Late_Pennsylvanian_300Ma.nc' ## Late Pennsylcanian ice, co2 200? 250?

1. input_dem="./indata/Map56_PALEOMAP_1deg_Early_Permian_295Ma.nc"

1. indatafile1='./indata/TOPicemsk.GLACD26kN9894GE90227A6005GGrBgic.nc'

1. input_dem="origodem.nc"
2. a_yr=14500

1. 1. load_glac1d_dem(indatafile1, input_dem, 14500)

1. 1. input one de scotese palaeomap dem!
2. def convert_to_t42(infilename1, outfilename1):

1. topo, lons, lats=convert_to_t21(input_dem, "demT21.nc")

1. topo, lons, lats=convert_to_t42(input_dem, "demT42.nc")

1. plt.imshow(topo,cmap='gist_earth')

1. plt.show()

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

1. a_resolution1="T42"

a_resolution1="T21" a_precision1=4 a_crashtolerant1=True a_landmap1="landmask.sra" a_topomap1="topo.sra"

1. 1. nowadays ca 0 BP
2. a_eccentricity1=0.01671022
3. a_obliquity1=23.44
4. a_lonvernaleq1=102.7
5. a_pCO21=360e-6

1. 1. 10000 yrs ago

1. a_eccentricity1=0.0194246086670259
2. a_obliquity1=24.230720588
3. a_lonvernaleq1=295.26651297
4. a_pCO21=265e-6

1. 1. 14500 yrs ago

1. a_eccentricity1=0.019595
2. a_obliquity1=23.6801
3. a_lonvernaleq1=221.5
4. (229.64+213.3)/2
5. a_pCO21=210e-6

1. 1. 25000 yrs ago

1. a_eccentricity1=0.0178681374211005
2. a_obliquity1= 22.408850897
3. a_lonvernaleq1=49.92
4. a_pCO21=180e-6

1. 1. cretaceous

a_eccentricity1=0.0167022 a_obliquity1=23.441 a_lonvernaleq1=282.7 a_pCO21=2500.0e-6

1. a_pCO21=700.0e-6

1. 1. early permian 295 ma
   2. late pennsylvanian 300 ma
2. a_eccentricity1=0.01671022
3. a_obliquity1=23.441
4. a_lonvernaleq1=282.7
5. a_pCO2=250.0e-6 ## ca 200 - 250 ppmvol
6. a_pCO21=180.0e-6
7. a_pCO21=100.0e-6

1. 1. permo-triassic boundary ca 250 ma
2. a_eccentricity1=0.01671022
3. a_obliquity1=23.441
4. a_lonvernaleq1=282.7
5. a_pCO21=1600.0e-6 ## cal1600 ppmvol 3000 ? 2000-4000

print("Exoplasim ...")

1. 1. if you run simu first time, you must set
2. firstrun=1

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(".")