1. "R" temperature downscaler

1. 1. cgenie muffin output && paleodem temperature downscaler

1. 1. 23.5.2023 v 0000.0004

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

library(pixmap) library(png) library(bmp) library(stringr)

library(OpenImageR) 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)

## invert #r=r*-1

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

#r[r[]>=1] <- #r[r[]<1] <- 1 #r[r[]==-2] <- 0

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)

}

year1=9

uxsize=720 uysize=360

1. 1. only input files

dataname1<-"./cgenie_1/fields_biogem_2d.nc" ## cgenie muffin output finename1<-"./maps/map445Ma.nc" ## paleodem file

1. 1. output files ...

coarsename1<-"atm_temp.nc"

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

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

edgename1="edge.png"

demname1= "planet_dem.nc"

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

1. edgename1="edges.png"

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, "atm_temp") # store the data in a 3-dimensional array dim(temp.array) nc_close(nc_data)

print(t)

temp.slice <- temp.array[, , year1]

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

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

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

plot(r)

writeRaster(r, "atm_temp.nc", "CDF", overwrite=TRUE)

1. stop(-1)

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

finer2<-finer1

finer2[finer2<0.0]<-0.0

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)

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

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

image(rout1)

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

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

finer3<-raster(finename1) finer4<-raster(finename1)

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)

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. *******************************************************************
2. GENIE COMPONENT SELECTION
3. *******************************************************************

ma_flag_ebatmos=.TRUE. ma_flag_goldsteinocean=.TRUE. ma_flag_goldsteinseaice=.TRUE. ma_flag_biogem=.TRUE. ma_flag_atchem=.TRUE.

1. 1. 2. cGENIE .config file parameter lines generated by muffingen v0.8 on: 191127 ###
2. INPUT FILE PATH

ea_1='./../../cgenie.muffin/genie-paleo/445eb05X' go_1='./../../cgenie.muffin/genie-paleo/445eb05X' gs_1='./../../cgenie.muffin/genie-paleo/445eb05X'

1. Grid resolution

GENIENXOPTS='$(DEFINE)GENIENX=36' GENIENYOPTS='$(DEFINE)GENIENY=36' GOLDSTEINNLONSOPTS='$(DEFINE)GOLDSTEINNLONS=36' GOLDSTEINNLATSOPTS='$(DEFINE)GOLDSTEINNLATS=36' GOLDSTEINNLEVSOPTS='$(DEFINE)GOLDSTEINNLEVS=16'

1. Topography

ma_fname_topo='445eb05X'

1. Assumed longitudinal offset of the grid

gm_par_grid_lon_offset=-180

1. Ocean depth scalar (m) [internally, parameter: dsc]

go_par_dsc=5000

1. Boundary conditions: EMBM

ea_topo='445eb05X' ea_taux_u='445eb05X.taux_u.dat' ea_tauy_u='445eb05X.tauy_u.dat' ea_taux_v='445eb05X.taux_v.dat' ea_tauy_v='445eb05X.tauy_v.dat' ea_adv_u='445eb05X.wvelx.dat' ea_adv_v='445eb05X.wvely.dat'

1. Boundary conditions: GOLDSTEIN

go_topo='445eb05X'

1. Boundary conditions: GOLDSTEIN sea-ice

gs_topo='445eb05X'

1. Boundary conditions: ALBEDO!

ea_par_albedo1d_name='445eb05X.albd.dat'

1. Boundary conditions: BIOGEM

bg_par_pindir_name='./../../cgenie.muffin/genie-paleo/445eb05X/' bg_par_windspeed_file='445eb05X.windspeed_uvaa.dat'

1. BIOGEM MISC
2. gas transfer coeff

bg_par_gastransfer_a=1.044

1. Grid resolution of solid Earth components

SEDGEMNLONSOPTS='$(DEFINE)SEDGEMNLONS=36' SEDGEMNLATSOPTS='$(DEFINE)SEDGEMNLATS=36' ROKGEMNLONSOPTS='$(DEFINE)ROKGEMNLONS=36' ROKGEMNLATSOPTS='$(DEFINE)ROKGEMNLATS=36'

1. Topography for solid Earth components

sg_par_pindir_name='./../../cgenie.muffin/genie-paleo/445eb05X/' sg_par_sed_topo_D_name='445eb05X.depth.dat' sg_par_sed_reef_mask_name='445eb05X.reefmask.dat' sg_par_sedcore_save_mask_name='445eb05X.sedcoremask.dat' rg_par_pindir_name='./../../cgenie.muffin/genie-paleo/445eb05X/' rg_topo='445eb05X.k1'

1. GEOLOGIC AGE DEPENDENT PARAMETERS
2. Solar constant (W m-2) ... don't forget to adjust it if not modern!!
   1. 1. ma_genie_solar_constant=1318.0
3. ... also, salinity should be set 1 PSU lower if it an ice-free World
   1. 1. go_saln0=33.9
4. Orbital parameters (modern, defaults)
   1. 1. ea_par_orbit_osce=0.0167 # eccentricity
      2. ea_par_orbit_oscsob=0.397789 # sine of obliquity
      3. ea_par_orbit_oscgam=102.92 # longitude of perihelion

ea_par_orbit_osce=0.06 # eccentricity ea_par_orbit_oscsob=0.374607 # sine of obliquity ea_par_orbit_oscgam=90 # longitude of perihelion

1. Ocean Mg/Ca concentrations (modern defaults, mol kg-1)

bg_ocn_init_35=1.028E-02 bg_ocn_init_50=5.282E-02

3. --- CLIMATE ---------------------------------------------------------
5. scaling for atmospheric CO2 radiative forcing, relative to 278 ppm

ea_radfor_scl_co2=0.5

2. --- BIOLOGICAL NEW PRODUCTION ---------------------------------------
4. No biology

bg_par_bio_prodopt='NONE'

2. --- DATA SAVING -----------------------------------------------------

bg_par_data_save_level=10 bg_ctrl_debug_lvl0=.true. ma_debug_loop=1

2. --- FORCINGS --------------------------------------------------------
4. add a geothermal heat flux (W m-2)
5. bg_ctrl_force_GOLDSTEInTS=.TRUE.
6. bg_par_Fgeothermal=100.0E-3
8. --- MISC ------------------------------------------------------------
11. --- END -------------------------------------------------------------