1. 3. Exoplasim planet running code
   4. exoplasim example
   5. stepper code
2. 2. 25.08.2023 0000.0003
3. 1. in ubuntu you must install
4. 1. pip3 install exoplasim[netCDF4]
   2. not
   3. "sudo pip3 install exoplasim[netCDF4]"

import numpy as np import math as math 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, seamasklevel1):

global NLAT global NLON

topo2=np.copy(topo)

seamasklevel2=seamasklevel1+1.0 topo2[topo2 < seamasklevel1] = seamasklevel1 masko=np.copy(topo2) masko[masko == seamasklevel1] = -9999999 masko[masko > seamasklevel1] = 1 masko[masko == -9999999 ] = 0

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, seamasklevel1):

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

dem000=nc['z'] dem=np.flipud(dem000) #dem=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) lts0=[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]

## lns0=[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]

lts1=np.array(lts0) lns1=np.array(lns0)

lns=lns1 lts=np.flip(lts1)

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) topo2=np.copy(topo) masko=np.copy(topo)

seamasklevel2=seamasklevel1+1.0

topo2[topo2 < seamasklevel1] = seamasklevel1

masko=np.copy(topo2) masko[masko == seamasklevel1] = -9999999 masko[masko > seamasklevel1] = 1 masko[masko == -9999999 ] = 0 #plt.imshow(demo) #plt.imshow(masko) #plt.imshow(topo2) #plt.show()

#grid=np.fliplr(masko) #def savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1): savenetcdf_single_frommem(outfilename1, "z", topo,lts,lns) savenetcdf_single_frommem("mapt21.nc", "z", topo2,lts,lns) savenetcdf_single_frommem("maskt21.nc", "z", masko,lts,lns) return(topo,lons,lats)

1. 1. stepper exoplasim ...

def run_exoplasim_wya(a_input_dem1,s_seamasklevel1, a_gridtype, a_layers, a_years,a_timestep,a_snapshots,a_ncpus,a_eccentricity,a_obliquity,a_lonvernaleq,a_pCO2,a_startemp,a_baseflux,a_yearlength,a_radius,a_gravity,a_rotationperiod):

print("Exoplasim runner ")

pressure1=1.0 a_pO2=(1-a_pCO2-0.79)*pressure1 a_pN2=(1-0.21-a_pCO2)*pressure1 a_pCO2=a_pCO2*pressure1

output_format=".nc"

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", a_seamasklevel1)

create_sras(topo, a_seamasklevel1)

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 } # !!!! fixed rotation! testplanet.configure( startemp=a_startemp, flux=a_baseflux,# Stellar parameters eccentricity=a_eccentricity, obliquity=a_obliquity, lonvernaleq=a_lonvernaleq, year=a_yearlength, fixedorbit=True, # Orbital parameters rotationperiod=a_rotationperiod, # Rotation IS FIXED!!! synchronous=False, topomap="topo.sra", landmap="landmask.sra", radius=a_radius, gravity=a_gravity, # 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()

looplen1=a_runsteps1

peen=0 runc1=1

print("Phase 2 !!! Stepper runner.")

for n in range(0,looplen1): print("Exoplasim runner year ",n)

a_years2=1 runc1=1

testplanet.run(years=1,crashifbroken=True) savename = 'planet_run_'+str(runc1) testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True) testplanet.save(savename) #tas=testplanet.inspect("tas") #tas2=np.ravel(tas) #tmean1=np.mean(tas2) #tmean2=math.floor( (tmean1-273.15)*100) #tmean3=tmean2/100 #print("Mean tas "+str(tmean3) ) tas=testplanet.inspect("tas") mint=testplanet.inspect("mint") maxt=testplanet.inspect("maxt") tas2=np.ravel(tas) tmean1=np.mean(tas2) tmean2=math.floor( (tmean1-273.15)*100) tmean3=tmean2/100 maxt2=np.ravel(maxt) maxt1=np.mean(maxt2) maxt2=math.floor( (maxt1-273.15)*100) maxt3=maxt2/100 mint2=np.ravel(mint) mint1=np.mean(mint2) mint2=math.floor( (mint1-273.15)*100) mint3=mint2/100 deltat3=maxt3-mint3 print("Mean tas "+str(tmean3) ) print("Max tas "+str(maxt3) ) print("Min tas "+str(mint3) ) print("Delta tas "+str(deltat3) )

print("Return.")

return(0)

print(" Exoplasim simulation code ---")

1. !!! snowball earth

1. input_dem="dema1.nc" ##dem of exoplanet
2. input_dem="elev1.nc" ##dem of exoplanet
3. input_dem="Map88_PALEOMAP_1deg_Cambrian_Precambrian boundary_540Ma.nc"

input_dem="Map01_PALEOMAP_1deg_Holocene_0Ma.nc"

a_modelname1="planet" a_workdir1="planet_run"

a_runsteps1=300 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" a_seamasklevel1=0.0

a_startemp=5772

1. a_baseflux=1365.2*0.92 # 650 ma

a_baseflux=1365.2 a_yearlength=365 a_radius=1.0 a_gravity=9.81

1. 1. 1. !!! ROTATION !!!!

a_rotationperiod=1

1. 1. earth nowadays

a_eccentricity1=0.016

a_obliquity1=23.44

a_lonvernaleq1=282

a_pCO21=1e-6

print("Exoplasim ...")

1. 1. attempt to run exoplasim stepper code

run_exoplasim_wya(input_dem, a_seamasklevel1, a_resolution1, a_layers1, a_years1,a_timestep1,a_snapshots1,a_ncpus1,a_eccentricity1,a_obliquity1,a_lonvernaleq1,a_pCO21,a_startemp,a_baseflux,a_yearlength,a_radius,a_gravity,a_rotationperiod)

print(".")