File:Gliese 12b example simulation sea land 1 r 1 1 1 1.png

https://github.com/alphaparrot/ExoPlaSim

Planet map generator in "diku" mapm 55555, adjusted curves with image processor: mone lower areas for sea

https://topps.diku.dk/torbenm/maps.msp

2. 1. exoplasim planet runner
3. 1. python 3, exoplasim
5. basic run settings T21
6. 1. 29.05.2024 0000.0002a3

import math import numpy as np import exoplasim as exo

timestep=30.0 years=100

Sol=1361.5

1. 1. basic input params

startemp=3296 # assumption Earth-like

1. 1. Gliese 12 b

planetname="Gliese 12 b" S1=1.63 flux=S1*Sol year=12.76140 rotationperiod=year

1. radius=0.958 ## rel to earth

radius=1.000 ## rel to earth

mass=math.pow(radius, 3.7)

eccentricity=0 obliquity=0

fixedorbit=True synchronous=True substellarlon=0

1. substellarlon=180

aquaplanet=False desertplanet=False vegetation=2 ## 0 no 1 static 2 dynamic stormclim=False aerosol=False co2weathering=False outgassing=0 evolveco2=False maxsnow=-1 ## False ? ozone=False vegaccel=1 seaice=False wetsoil=False

glaciers1= { "toggle": True, "mindepth":2, "initialh":-1 }

1. landmap="Alderaan_surf_0129.sra"
2. topomap="Alderaan_surf_0172.sra"

1. landmap=None
2. topomap=None

landmap="mapa_surf_0129.sra" topomap="mapa_surf_0172.sra"

1. relCO2=0.95
2. relN2=1-relCO2

relCO2=280e-6 relO2=0.21 relN2=1-relCO2-relO2

1. relCO2=0.95
2. relO2=0.05
3. relN2=1-relCO2-relO2

density=mass/np.power(radius,3)*5.519 vesc=np.sqrt(mass/radius)*11.186 geese=mass/(radius*radius) geeseg=geese*9.80665

patm1=geese*mass*radius*radius ## estimation patm2=math.pow(radius, 2.4)

gravity=geeseg

1. pressure=(patm1+patm2)/2
2. pressure=10
3. pressure=0.1

pressure=1

pN2=pressure*relN2 pCO2=pressure*relCO2 pO2=pressure*relO2

print(" Name ", planetname) print(" insol S0 ",round(S1,3) ) print(" insol Wm2 ",round(flux,2) ) print(" mass ",round(mass,2) ) print(" radius ",round(radius,2) ) print(" radius km ",round((radius*6371),2) ) print(" density ",round(density,2) ) print(" vesc ",round(vesc,2) ) print(" gee_e ",round(geese,2) ) print(" gee_ms2 ",round((geese*9.81),2) ) print(" P_atm ",round(pressure,6) ) print(" O2: ",relO2 ,"N2: ",relN2 ,"CO2: ",relCO2);

print("Running wait long time ...")

1. quit(-1)

planeta = exo.Model(workdir="planeta_run",modelname="Planeta",ncpus=4,resolution="T21",outputtype=".nc")

planeta.configure(year=year, wetsoil=wetsoil,pO2=pO2, ozone=ozone, vegaccel=vegaccel, seaice=seaice, glaciers=glaciers1, maxsnow=maxsnow, evolveco2=evolveco2, outgassing=outgassing, co2weathering=co2weathering, vegetation=vegetation, stormclim=stormclim,landmap=landmap, topomap=topomap,startemp=startemp, flux=flux, eccentricity=eccentricity,obliquity=obliquity,fixedorbit=fixedorbit,synchronous=synchronous,substellarlon=substellarlon,rotationperiod=rotationperiod,radius=radius,gravity=gravity,aquaplanet=aquaplanet,desertplanet=desertplanet,pN2=pN2,pCO2=pCO2,timestep=timestep,snapshots=False,physicsfilter="gp|exp|sp")

print(" Run, wait a long time...") planeta.run(years=years,crashifbroken=True)

planeta.exportcfg()

planeta.run(years=years,crashifbroken=True)

planeta.finalize("Planeta",allyears=True,keeprestarts=True)

planeta.save()

1. 1. diku "planet map generator" map converter to exoplasim input test
2. maps from https://topps.diku.dk/torbenm/maps.msp
   1. must use here grayscale palette outrput option !

1. 1. normalizes map to 0...1, then optional filtering, then normalize to depest and highest points
   2. python3 source code
2. 1. 29.5.2024 v 0000.0001b

import os import math as math import numpy as np

from scipy.interpolate import interp2d

import matplotlib.pyplot as plt

1. from skimage.transform import resize

import skimage as ski

from PIL import Image, ImageFilter

from imageio import imwrite import imageio

import netCDF4

1. import exoplasim as exo

1. 1. main ctrls

top1=7000 floor1=-8000

sealevel1=0.08 ## must use here small value! level in original bittmap 0...1 iname1="Map-555.bmp"

imagetype1=1 ## rgb

filtermap1=1 ## 1: np.power(map1a,math.exp(1)*2 above sealevel

NLAT=0 NLON=0

def convert_image_to_exoplasim_maps_1(iname1, sealevel1, top1, floor1,imagetype1, filtermap1):

imin1=Image.open(iname1) img1 = np.asarray(imin1).astype(float) max1=np.max(img1) min1=np.min(img1) del1=max1-min1 map0=(img1-min1)/del1

map1a=np.copy(map0) #map1=np.exp(map0)/np.exp(1) #map1a=np.power(map1a, 15) if(filtermap1==1): #map1a=np.where(map1a>sealevel1, np.exp(np.exp(map1a)), map1a) map1a=np.where(map1a>sealevel1, np.power(map1a,math.exp(1)*2), map1a)

max2=np.max(map1a) min2=np.min(map1a) del2=max2-min2 map1=(map1a-min2)/del2

if(imagetype1==1): map1=map1[:,:,0] else: map1=map1

shape1=np.shape(map1) width1=shape1[1] height1=shape1[0] imsize1=width1*height1 #quit(-1) map2=np.copy(map1) map2=np.where(map2<sealevel1,0,map2) map2=np.where(map2<sealevel1,map2*floor1,map2*top1)

landpixels1=np.count_nonzero(map2) landpercent1=round( ((landpixels1*100)/imsize1),1) landmean1=np.mean(map2)

print(width1, height1) print(" Land % " , landpercent1) print(" Land mean height m:", landmean1) #map2a=(np.abs(map2)*256*256).astype(np.uint16) map2a=map2 #print (np.max(map2a)) #quit(-1) #print(map2a) #quit(-1) imout1 = Image.fromarray(map2a) imageio.imwrite('map2.png', map2a.astype(np.uint16)) #imout1.save("map1.png") #quit(-1) map3=np.roll(map2a, int(width1/2), axis=1)

#plt.imshow(map3) #plt.show() #imout2 = Image.fromarray(map3) imageio.imwrite('map1.png', map3.astype(np.uint16)) mask2=np.copy(map2a) mask3=np.copy(map3) mask2=np.where(mask2<1,0, 65535) mask3=np.where(mask3<1,0, 65535) #imout3 = Image.fromarray(mask2) #imout4 = Image.fromarray(mask3) imageio.imwrite('mask2.png', mask2.astype(np.uint16)) imageio.imwrite('mask1.png', mask3.astype(np.uint16)) image1 = Image.open("mask1.png") image1 = image1.convert("L") image1 = image1.filter(ImageFilter.FIND_EDGES) image1.save("edge1.png") image2 = Image.open("mask2.png") image2 = image2.convert("L") image2 = image2.filter(ImageFilter.FIND_EDGES) image2.save("edge2.png") savenetcdf2("map1.nc", "z", map2a) conv_t21("map1.nc", "inmap", 0) return(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 conv_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) indimx=len(inlons) indimy=len(inlats)

latlen=len(inlats) lonlen=len(inlons)

#print(lonlen, latlen)

indimx=lonlen indimy=latlen

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

demo=ski.transform.resize(dem2, (shapey, shapex))

#topo=np.flipud(demo) topo=np.copy(demo) topo2=np.copy(topo) masko=np.copy(topo)

#plt.imshow(topo) #plt.show() #quit(-1)

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)

writeSRA("mapa",129,np.flipud(topo2),NLAT,NLON) writeSRA("mapa",172,np.flipud(masko),NLAT,NLON) return(topo,lons,lats)

def savenetcdf2(outfilename1, outvarname1, map1):

shape1=np.shape(map1) width1=shape1[1] height1=shape1[0] xoutlons1=np.linspace(0,360,width1) xoutlats1=np.linspace(-90,90, height1)

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[:, :] =np.flipud(map1[:]) ncout1.close() return 0

1. 1. main program

convert_image_to_exoplasim_maps_1(iname1, sealevel1, top1, floor1, imagetype1, filtermap1)

print(".")

quit(-1)