File:Teegarden b temperature if locked earthlike 3mearths patm 9 o2 21 1.png

1. https://www.aanda.org/articles/aa/full_html/2024/04/aa48033-23/aa48033-23.html

1. https://www.aanda.org/articles/aa/pdf/2024/04/aa48033-23.pdf
2. A&A, 684, A117 (2024)
3. Teegarden’s Star revisited
4. A nearby planetary system with at least three planets★,★★
5. S. Dreizler tet al

2. 1. exoplasim planet runner
   2. python 3, exoplasim
4. basic run settings T21
5. 1. 29.10.2023 0000.0002a2

import math import numpy as np import exoplasim as exo

timestep=30.0 years=100

Sol=1361.5

1. 1. basic input params

1. https://www.aanda.org/articles/aa/full_html/2024/04/aa48033-23/aa48033-23.html

1. https://www.aanda.org/articles/aa/pdf/2024/04/aa48033-23.pdf
2. A&A, 684, A117 (2024)
3. Teegarden’s Star revisited
4. A nearby planetary system with at least three planets★,★★
5. S. Dreizler tet al

startemp=3034 ## Teegarden's Star Zec 2024

1. 1. mstar 0.097
   2. rstar 0.120
   3. feh -0.11
   4. star prot 96.2

planetname="Teegarden b"

1. S1=7.22e-6/math.pow(0.0259,2)
   1. tew 277+-5

S1=1.08 flux=S1*Sol year=4.90634 rotationperiod=year

1. mass=1.16*1.0

mass=3 radius=math.pow(mass, 0.27)

1. radius=1.04
2. mass=math.pow(radius, 3.7)

eccentricity=0.03 obliquity=0

fixedorbit=True synchronous=True

1. 1. substellarlon=0

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=True wetsoil=True

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.96
2. relN2=1-relCO2
3. relO2=0.00

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

1. relCO2=0.95
2. relO2=0.045
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

pressure=math.pow(mass,2)*1 ## often assumed m2

1. pressure=math.pow(mass,2.3) ## some boosted effect ...
2. pressure=math.pow((mass),2)*6e-3 ## thin Mars-like atmos!

1. pressure=10
2. pressure=0.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(" Patm ",round(pressure,6) )

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

1. quit(-1)

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

1. 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, #aerosol=aerosol,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")

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. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3
2. 1. diku "planet map generator" map converter to exoplasim input test
3. maps from https://topps.diku.dk/torbenm/maps.msp
   1. must use here grayscale palette output option !

1. 1. normalizes map to 0...1, then np.exp(z), then normalize to depest and highest points
   2. python3 source code
2. 1. 24.6.2024 v 0000.0000c2

import os import math as math import numpy as np

from scipy.interpolate import interp2d import matplotlib.pyplot as plt

from PIL import Image, ImageFilter

from imageio import imwrite import imageio

import netCDF4

1. import exoplasim as exo

1. 1. main ctrls

iname1="Map-777.bmp"

top1=10000 floor1=-10000

land_height_scaling_exponent_1=2 ## land height scaling exponent 6

sealevel1=0.292 ## if land scaling exponent hight mush use here small value! level in original bittmap 0...1

def t21_meanz(inmap):

   latitudes0=np.array([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])

   radlats1=np.radians(latitudes0)
   coslats1=np.cos(radlats1)
   height1=len(radlats1)
   seas1=0
   all1=0
   landmeanz1=0

   for iy in range(0,(height1)):
       line1=inmap[iy,:]
       lek1=abs(coslats1[iy])
       seadd0=np.count_nonzero(line1 == 0)
       #print(seadd0)
       seadd1=seadd0/len(line1)
       seadd2=lek1*seadd1
       seas1=seas1+seadd2
       all1=all1+lek1
       lamez0=line1[line1>0].mean()
       if (math.isnan(lamez0)): lamez0=0
       lamez1=lamez0*lek1
       landmeanz1=landmeanz1+lamez1

   landmeanz2=round((landmeanz1/height1),2)

   seapercent2=round(((seas1*100)/all1),2)
   landpercent2=100-seapercent2
   #print(seas1)
   #print(all1)
   #print(seapercent2)
   #print(landpercent2)
   return(landpercent2, seapercent2, landmeanz2)

def count_spherical_land_sea_meanz(inmap):

   radlats1=np.linspace(-np.pi/2,np.pi, height1)
   coslats1=np.cos(radlats1)
   seas1=0
   all1=0
   landmeanz1=0

   for iy in range(0,(height1)):
       line1=inmap[iy,:]
       lek1=abs(coslats1[iy])
       seadd0=np.count_nonzero(line1 == 0)
       #print(seadd0)
       seadd1=seadd0/len(line1)
       seadd2=lek1*seadd1
       seas1=seas1+seadd2
       all1=all1+lek1
       lamez0=line1[line1>0].mean()
       if (math.isnan(lamez0)): lamez0=0
       lamez1=lamez0*lek1
       landmeanz1=landmeanz1+lamez1

   landmeanz2=round((landmeanz1/height1),2)

   seapercent2=round(((seas1*100)/all1),2)
   landpercent2=100-seapercent2
   #print(seas1)
   #print(all1)
   #print(seapercent2)
   #print(landpercent2)
   return(landpercent2, seapercent2, landmeanz2)

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) 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) landpercent2, seapercent2, landmeanz2=t21_meanz(topo2) print("") print("Output T21 dem \n Land % ", landpercent2, "\n Mean land z ", landmeanz2) #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

NLAT=0 NLON=0

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

1. map1=np.exp(map0)/np.exp(1)

map1a=np.copy(map0)

1. 1. land_height_scaling_exponent_1=15

map1a=np.power(map1a, land_height_scaling_exponent_1)

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

map1=(map1a-min2)/del2

map1=map1[:,:,0]

shape1=np.shape(map1)

width1=shape1[1] height1=shape1[0]

imsize1=width1*height1

1. quit(-1)

map2=np.copy(map1) map3=np.copy(map1)

map2=np.where(map2<sealevel1,0,map2)

map2=np.where(map2<sealevel1,map2*floor1,map2*top1)

1. 1. abyss mat too ...

map3=np.where(map3<sealevel1,map3*floor1,map3*top1)

landpixels1=np.count_nonzero(map2)

landpercent1=round( ((landpixels1*100)/imsize1),1)

1. quit(-1)

landpercent2, seapercent2, landmeanz2=count_spherical_land_sea_meanz(map2)

1. map2a=(np.abs(map2)*256*256).astype(np.uint16)

map2a=map2

1. print (np.max(map2a))

1. quit(-1)

1. print(map2a)

1. quit(-1)

imout1 = Image.fromarray(map2a)

imageio.imwrite('map2.png', map2a.astype(np.uint16))

1. imout1.save("map1.png")

1. quit(-1)

map3=np.roll(map2a, int(width1/2), axis=1)

plt.imshow(map3) plt.show()

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

1. imout3 = Image.fromarray(mask2)
2. 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) print() print() print() print("Output coarse map params:") conv_t21("map1.nc", "inmap", 0)

1. print(width1, height1)
2. print(" Land %" , landpercent1)

print() print("Accurate dem") print(" Land  % ", landpercent2) print(" Sea  % ", seapercent2)

print(" Land mean height ", landmeanz2)

print(".")

1. "R" temperature downscaler

1. 1. exoplasim output && paleodem temperature 6.5 C/1km downscaler

1. 1. 4.6.2023 v 0000.0001a

library(stringr)

library(raster) library(ncdf4) library(rgdal) library(viridis)

library(pixmap) library(png) library(bmp) library(OpenImageR)

1. library(graphx)

save_as_nc<-function(outname1,r, outvar1, longvar1, unit1) { ext2<-c(-180,180, -90,90) extent(r)<-ext2 crs(r)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(r, outname1, overwrite=TRUE, format="CDF", varname=outvar1, varunit=unit1,

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

}

create_mask_rasters<-function(infile1,maskname1,demname1) {

ur1<-raster(infile1) ur1[ur1[]<1] <- 0 ## jn warning flip! #rdem1=flip(ur1) r=flip(ur1) #r=ur1 rdem1=ur1

mini=minValue(rdem1) maxi=maxValue(rdem1) delta=maxi-mini rdem2=(rdem1/delta)*255

dims<-dim(r)

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

zeros1<-sum(r[]==0) ones1<-sum(r[]==1)

print(zeros1) print(ones1) #stop(-1)

image(r)

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

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

mask0<-r mask1<-mask0[]

idem1<-rdem2[] mask2<-matrix(mask1, ncol=cols, nrow=rows ) idem2<-matrix(idem1, ncol=cols, nrow=rows ) #mask3<-t(mask2) #idem3<-t(idem2) mask3<-t(mask2) idem3<-t(idem2) mask4 <- apply(mask3, 2, rev) idem4 <- apply(idem3, 2, rev) r <- writePNG(mask4, maskname1) r <- writePNG(idem4, demname1) }

createmask <- function(dem1, demname2, maskname1, maskname2,maskname3, maskname4, masksealevel1) {

 dem2<-(dem1+abs(masksealevel1))
 
 mask2<-dem2
 mask1<-dem2
 
 mask2[mask2 > 1 ]=255
 mask2[mask2 < 0 ]=0 
 
 mask1[mask1 > -1 ]=0
 mask1[mask1 < 0 ]=255  
 
 	zeros1<-sum(mask2[]==0)

ones1<-sum(mask2[]==255) all1=zeros1+ones1

landsearatio1=(100*zeros1)/all1

print(zeros1) print(ones1) print (all1) print(800*1600) print("Sea per cent of total surface:") print(landsearatio1)

#print("TEBUK STOP !!!") #stop(-1)

 image(dem1)
 image(mask1)
 image(mask2)
   
 #dem2<-flip(dem2, direction= "y")
 #mask1<-flip(mask1, direction="y") 
 #mask2<-flip(mask2, direction="y")    
 
 ext1<- extent(0, 360, -90, 90)
 ext2<- extent(-180, 180, -90, 90) 
 
 mask30<-mask1
 extent(mask30) <- ext1
 mask3 <- rotate(mask30)
 extent(mask3) <- ext1
 mask40<-mask2
 extent(mask40) <- ext1
 mask4 <- rotate(mask40)
 extent(mask4) <- ext1
 

1. save_as_nc(demname2,dem2, "Band1", "Band1", "")
2. save_as_nc(maskname1,mask1, "Band1", "Band1", "")
3. save_as_nc(maskname2,mask2, "Band1", "Band1", "")
4. save_as_nc(maskname3,mask3, "Band1", "Band1", "")
5. save_as_nc(maskname4,mask4, "Band1", "Band1", "")

}

downscale_temperature_65 <- function (coarse_rastera, fine_rastera) { ## methods: 0 delta, 1 spatialeco, 2 dissever, 3 temperature lapse 6.5 C/1 km lm

   print ("Downscaling temperature, 6.5 C/1 km")			

coarse_raster<-coarse_rastera fine_raster<-fine_rastera p1<-fine_raster p2<-fine_raster

plot(fine_raster) plot(coarse_raster, col=viridis(200))

coarseoro<- resample(p1, coarse_raster) coarseoro_big<-resample(coarseoro, p1) orodelta<-(p1-coarseoro_big)

baset1 <- resample(coarse_raster, p1)

raster_stack <- stack(p1,p2)

orotemp<-orodelta*0.0065*-1

coarseorotemp<- resample(orotemp, coarse_raster) coarseorotemp_big<-resample(coarseorotemp, p1)

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

#outtempr<-rotate(outtemp)

#plot(outtempr)

     return(outtemp)

}

1. 1. warning: check, if you need
   2. flip input MOST raster, You can do this witk "lsm" raster

1. 1. only input files
2. dataname1<-"../tea_locked_earthlike_mca1p16_pca1p3_co295_1/MOST.00026.nc"

dataname1<-"../lockead_easuperearthlike_mass_3me_1/MOST.00027.nc"

finename1<-"../map1.nc"

1. 1. output files ...

coarsename1<-"tassi.nc"

maskname1= "mask_land.png" maskname2= "mask_sea.png"

edgename1="edge.png"

1. maskname3= "mask_land_rotated.nc"
2. maskname4= "mask_sea_rotated.nc"

demname1= "planet_dem.nc"

uxsize=720 uysize=360

ext1<-extent(-180,180,-90,90) crs1<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

1. stop(-1)

nc_data <- nc_open(dataname1)

lon <- ncvar_get(nc_data, "lon") lat <- ncvar_get(nc_data, "lat", verbose = F) t <- ncvar_get(nc_data, "time")

temp.array <- ncvar_get(nc_data, "tas") # store the data in a 3-dimensional array mask.array <- ncvar_get(nc_data, "lsm") ## for ratster flip check only

1. plot(mask.array)

1. stop(-1)

temp.array <-temp.array-273.15

dim(temp.array) nc_close(nc_data)

print(t)

1. temp.slice <- temp.array[, , 6]

print (dim(temp.array))

temp.slice<-apply(temp.array, c(1,2), mean) mask.slice<-apply(mask.array, c(1,2), mean)

1. plot(temp.slice)

1. stop(-1)

r <- raster(t(temp.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

1. plot(r)

maskr <- raster(t(mask.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

plot(maskr)

1. stop(-1)

1. 1. check whethet must fli ras

1. 1. r <- flip(r, direction='y')

extent(r)<-ext1 crs(r)<-crs1 extent(maskr)<-ext1 crs(maskr)<-crs1

1. plot(r)

1. plot(maskr)

writeRaster(r, coarsename1, "CDF", overwrite=TRUE)

1. stop(-1)

coarser1<-raster(coarsename1) finer1<-raster(finename1)

extent(finer1)<-ext1

plot(coarser1)

1. plot(finer1)

1. stop(-1)

finer2<-finer1

finer2[finer2<0.0]<-0.0

plot(coarser1) plot(finer2)

1. stop(-1)

matrix1 <- matrix(1:(uxsize*uysize), nrow=uysize, ncol=uxsize)

sabluna1 <- raster(matrix1) extent(sabluna1) <- ext1 crs(sabluna1)<-crs1

finer4<-resample(finer2, sabluna1) finer0<-finer4

rout1<-downscale_temperature_65(coarser1, finer4)

plot(rout1)

1. stop(-1)

crs(rout1)<-crs1 extent(rout1 )<-ext1

1. image(rout1)

writeRaster(rout1, "temp_dskaled1.nc", overwrite=TRUE, format="CDF", varname="atm_temp", varunit="degC",

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

finer3<-finer0

finer3[finer3>0]<-255 finer3[finer3<=0]<-0

writeRaster(finer3, "mask1.tiff", overwrite=TRUE, format="GTiff", varname="mask", varunit="unit",

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

dem1<-raster(finename1)

writeRaster(dem1, "demo1.nc", overwrite=TRUE, format="CDF", varname="mask", varunit="unit",

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

masksealevel1=0

create_mask_rasters(finename1,maskname1,demname1)

image1 = readImage(maskname1)

1. image1

1. print (dim(image1))

image2=image1

1. str(image2)

1. stop(-1)

res1 = edge_detection(image2, method = 'Frei_chen', conv_mode = 'same')

1. plot(res1)

res2=res1

res2[res2[]==0]=2 res2=res2-1

writeImage(res2, edgename1)

image2[image2[]==0]=2 image2=image2-1

writeImage(image2, maskname2)

1. createmask(dem1,demname1,maskname1, maskname2, maskname3, maskname4, masksealevel1)