File:Trias 240ma 1.png

1. 1. process dem file to mask
   2. and flatten sea

install_packages=0

if(install_packages==1) {

install.packages("raster") install.packages("ncdf4") install.packages("rgdal") install.packages("png") }

library(raster) library(ncdf4) library(rgdal) library(png)

file1="./trias240_60.nc" file2="dem.nc" file3="dem.tif" outmask1="out_mask.png"

ur1<-raster(file1)

ur1[ur1[]<1] <- 0

1. image(ur1)

1. plot(ur1)

lonr1 <- init(ur1, 'x') latr1 <- init(ur1, 'y')

crs(ur1)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

writeRaster(ur1, file2, overwrite=TRUE, format="CDF", varname="Band1", varunit="m",

       longname="Band1", xname="lon",   yname="lat")

writeRaster(ur1, file3, overwrite=TRUE, format="GTiff", varname="Band1", varunit="m",

       longname="Band1", xname="lon",   yname="lat")

crs(lonr1)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" crs(latr1)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

writeRaster(lonr1, "lons.nc", overwrite=TRUE, format="CDF", varname="Band1", varunit="deg",

       longname="Band1", xname="lon",   yname="lat")
       

writeRaster(latr1, "lats.nc", overwrite=TRUE, format="CDF", varname="Band1", varunit="deg",

       longname="Band1", xname="lon",   yname="lat")
              
     
     
     
       

r=ur1

dims<-dim(r)

dims

r[r[]<1] <- 0 r[r[]>0] <- 1

image(r)

1. stop(-1)

print (dims[1]) print (dims[2])

rows=dims[2] cols=dims[1]

1. stop(-1)

mask0<-r

mask1<-mask0[]

mask2<-matrix(mask1, ncol=cols, nrow=rows )

mask3<-t(mask2)

r <- writePNG(mask3, outmask1)

plot(r)

1. png('mask.png', height=nrow(r), width=ncol(r))
   1. plot(r, maxpixels=ncell(r))
2. image(r, axes = FALSE, labels=FALSE)
3. dev.off()

1. 3. Ecoplasim planet running code
   4. exoplasim example
   5. in ra
   6. convert to T21, input netcdf
   7. load one lon, lat, z grid
   8. or Tarasov glac1d grid
   10. 08.01.2022 0000.0005
2. 1. MPI NOTE: if you use more than
   2. one processor, you cannot in most cases run MPI in root
   3. in ubuntu you must install
3. 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 run_exoplasim_b(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=".npz" 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.Model(workdir="testplanet_run",modelname="TESTPLANET",ncpus=a_ncpus,resolution="T21") #testplanet= exo.Earthlike(workdir="planet_run",modelname="PLANET",ncpus=a_ncpus,resolution="T21",outputtype=output_format, crashtolerant=True) testplanet= exo.Earthlike(workdir="planet_run",modelname="PLANET",ncpus=a_ncpus,resolution=a_gridtype,layers=a_layers, outputtype=output_format, crashtolerant=True)

## earth 21000 BP glaciers1= { "toggle": True, "mindepth":2, "initialh":-1 }

testplanet.configure( startemp=5772.0, flux=1367,# 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, # Bulk properties #seaice=False, #maxsnow=False, #glaciers=False, #stormclim=False, #vegetation=0, wetsoil=True, #alters albedo of soil based on how wet it is

               vegetation=2,                               #toggles vegetation module; 1 for static vegetation, 2 to allow growth
               vegaccel=1, 

seaice=True, maxsnow=-1, glaciers=glaciers1, #stormclim=True, #vegetation=0, 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") # Model dynamics

testplanet.exportcfg()

print("Running ExoPlasim ... ")

print("K-Pg limit transient test.")

#testplanet.run(years=a_years,crashifbroken=True)

a_years=50 runc1=1

print("Phase 1 !!! Init run. Summer.") testplanet.run(years=a_years,crashifbroken=True)

savename = 'planet_run_'+str(runc1) testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True) testplanet.save(savename)

print("Phase 2 !!! Impact, Winter, springwinter ...")

fluxes1=[1367*1.0,1367*1.0] looplen1=len(fluxes1)

peen=0 runc1=2

print ("Impact flux", fluxes1[0])

for n in range(0,looplen1): darkflux=fluxes1[peen] testplanet.modify(flux=darkflux) #number of output times (months) in the output files testplanet.exportcfg() a_years2=1 runc1=2 print("Winter: Year from impact ",n) testplanet.run(years=1,crashifbroken=True) savename = 'planet_run_'+str(runc1) testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True) testplanet.save(savename)

print("Phase 2 !!! Winter, Spring ...") darkflux=1367 testplanet.modify(flux=darkflux) #number of output times (months) in the output files testplanet.exportcfg() a_years2=50 runc1=3 print("Spring ...") testplanet.run(years=a_years,crashifbroken=True)

savename = 'planet_run_'+str(runc1) testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True) testplanet.save(savename)

print("Return.") return(0)

print(" Exoplasim simulation code ---")

1. 1. jn warning maybe nok

1. input_dem='./indata/indem.nc'

input_dem="./indata/Map16_PALEOMAP_1deg_KT_Boundary_65Ma.nc"

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

input_dem="./Scotese_Wright_2018_Maps_1-88_1degX1deg_PaleoDEMS_nc_v2/Map47_PALEOMAP_1deg_Middle_Triassic_240Ma.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="trias" a_workdir1="trias_run"

a_runsteps1=51 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. triassic middle

a_eccentricity1=0.0167022 a_obliquity1=23.441 a_lonvernaleq1=102.7

1. a_pCO21=900.0e-6
2. a_pCO21=500.0e-6
3. a_pCO21=1200.0e-6

a_pCO21=1400.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=102.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=102.7
5. a_pCO21=1600.0e-6 ## cal1600 ppmvol 3000 ? 2000-4000

print("Exoplasim ...")

run_exoplasim_b(input_dem, a_resolution1, a_layers1, a_years1,a_timestep1,a_snapshots1,a_ncpus1,a_eccentricity1,a_obliquity1,a_lonvernaleq1,a_pCO21)

print(".")