1. 3. Hadcm3b 60kyr netcdf downscaler
   5. Peltier Glac1D Etopo1
2. 1. Language "R" 4.0.2
   3. v 1015.05
   4. 28.9.2020
3. 1. Externel requirements:
4. You need in most cases Etopo1, Peltier and Tarasov Glac 1d dataset to
5. run this script
6. 1. WARNING1: Maybe inaccuracy due to inaccurate input orography data
   4. WARNING2: plotting is under development stage, must visualize w/Panoply or other tools

install_libraries=FALSE

if(install_libraries==TRUE) {

install.packages("raster")
install.packages("rgdal")
install.packages("sp")
install.packages("spatialEco")
install.packages("ncdf4")
install.packages("dissever")
install.packages("viridis")
install.packages("dplyr")
install.packages("lattice")
install.packages("RColorBrewer")
install.packages("pals")
install.packages("rgeos")
install.packages("sp")
install.packages("reshape2")
install.packages("data.table")
install.packages("stringr")
install.packages("rlist")
install.packages("pipeR")
install.packages("maptools")
install.packages("gdata", dependencies=TRUE)
install.packages("abind")
install.packages("Cairo")
install.packages("Rfast")

}

library(raster) library(rgdal) library(ncdf4) library(lattice) library(maptools) library(rgeos) library(spatialEco) library(dissever)

library(RColorBrewer) library(viridis) library(pals) library(data.table) library(stringr) library(rlist) library(pipeR) library(Rfast)

library(sp) library(reshape2)

library(dplyr)

1. library(gdata)

library(abind)

library(Cairo)

1. 1. NOTE SET THESE

1. 1. set orografy creating off

1. argreading=0 ## WARNING NOT IMPLEMENTED default 0
2. kalk_distance=0 # for speed 0, if u not neet to calc distance
3. kalk_tables=0 # suppress erroris from tables, id kalk_distance=0

1. 1. must be 1, if you create output foor certain year, first

1. 1. Control variables

arg_reading=1

orography_preprocess=1

1. 1. get dordogne srtm

get_srtm_data=0

1. 1. gtopo30

get_gtopo30_data=0

1. 2. output area type
   3. 0: normal local output
   4. 1: global output, according to Glac!D
   5. NOT OK 2: custom orography, local NOK

1. 1. acquire trace21ka temperature data, default 1

capture_temperature=0 capture_rainfall=1

1. 1. na fill method 0,1, 2 or 4
   2. method 4 resample
   3. method 2 is very slow

hadcm3_60ka_seafill_method=1

1. 1. if method 4, sample cols and rows
2. samplecols1=50
3. samplerows1=25

samplecols1=75 samplerows1=35

1. 1. acquire hadcm3 60 kyr data
   2. Warning: you must edit procedure to acquire right dataser

capture_hadcm3_60ka_data=1

1. 1. acquire trace 26k data

capture_trace26k_data=0

1. 1. acquire glac1d elevation data,
   2. default 2
   3. 1 Tarasov glac1d (26 ka)
   4. 2 Peltier 1ce 6G-D (122 ka) NOTE: must use, when age over 26ka!

capture_elevation=2

1. 1. generate rain shadows, assume westerly wind by default

make_rain_shadows=1

calculate_distance=0

1. 1. either global_ds or normal_ds must set 1
   2. normal etopo1 area downscaling normal_ds=3
2. 1. 1: temperature downscaling: normal_ds=1
   2. 3: rainfall downscaling: normal_ds=3

normal_ds=3

1. 1. either global_ds or normal_ds must set 1
   2. global area downscaling to glac1d

global_ds=0

1. ! afterburner accurate "kustomorofilee1.nc" downscaling, default 0
   1. need accurate orog file (GTOPO30 slice, SRTM etc.)

accurate_ds=0

1. kustomorofilee1="./predata/europe30.nc"
   1. plotting is under alpha stage!

plot_result=0

1. 1. downscaling method from orography
   2. 0 delta, 1 spatialeco, 2 dissever, 3 temperature delta lapse rate 6,5 c/ 1 km
   3. method 2 slow

method_ds_oro=1

1. 1. downscaling method for temperature

method_ds_var=1

kustomorofilee1="./predata/dordogne.nc"

1. kustomorofilee1="./predata/europe30.nc"

1. 1. year to render

fage=40740

1. 1. fage=13000
2. fage=44950
3. fage=38000

1. 1. month of year to render
   2. month 0 sum of all monnths: tex annual rainfall
   3. month 13 average of all months

fmonth=0

1. 1. age in tarasov ice sheet

fage2=fage

2. fage=56000
   1. number of years to average, since "fage"

numyears=20

1. 1. sea level, below current: eq 120 means height of -120 to current
   2. not needed set by default, calculated from age w/ thisscript

sealevel=-70

out_text3="X" out_text4= "Europe, "

if(capture_temperature==1) { out_text_1="TS" out_text_2="°C"

if(fmonth==7) { out_text_3=paste0("Temperature of July, ",out_text4) }

if (fmonth==1) { out_text_3=paste0("Temperature of January, ",out_text4) } }

if(capture_rainfall==1) { out_text_1="Precip" out_text_2="mm"

if(fmonth==7) { out_text_3=paste0("Precipitation of July, ",out_text4) }

if (fmonth==1) { out_text_3=paste0("Precipitation of July, ",out_text4) }

if(fmonth==0) { out_text_3=paste0("Annual precipitation , ",out_text4) } }

out_text_3=paste0(out_text_3,toString(fage)) out_text_3=paste0(out_text_3, " years ago")

1. out_text_1="Precipitation"
2. out_text_2="mm"
3. out_text_3="Annual precipitation, Europe, "
4. out_text_3=paste0(out_text_3,toString(fage))
5. out_text_3=paste0(out_text_3, " years ago")

1. 1. area to render, lon1, lon2 lat1 lat2 rectangle

1. 1. europe

lon1=-15.0 lon2=40.0 lat1=30.0 lat2=70.0

1. beringia
2. lon1=140
3. lon2=240
4. lat1=50.0
5. lat2=80.0

1. east asia
2. lon1=80
3. lon2=200
4. lat1=40.0
5. lat2=80.0

1. 1. global
2. klopaali1=1

1. lon1=-180.0
2. lon2=180.0
3. lat1=-90.0
4. lat2=90.0

1. 1. dordogne 1

1. lon1=-2.0
2. lon2=4.0
3. lat1=42.0
4. lat2=48.0

1. 1. dordogne 2
2. lon1=0.0
3. lon2=2.0
4. lat1=44.0
5. lat2=45.5

1. 1. dordogne 3
2. lon1=-1.0
3. lon2=2.0
4. lat1=44.0
5. lat2=45.5

if(fage>25999) {

1. NOTE: must use, when age over 26ka!

if(capture_elevation==1) { capture_elevation=2 } }

if(normal_ds==1) { if(capture_rainfall==1) { normal_ds=3 make_rain_shadows=1 } }

filename1="KKK" tslocation1=-999 elevlocation1=-999

1. 1. calculate location in tarasov ice sheet data

elevlocation1<<-round( (fage2-26000)/100)*-1

variaabeli="TS" unitti="Celcius" lonkkunimi="Temperature" lonvar1="lon" latvar1="lat"

1. 1. hadcm3b 60ka files path
2. hadbasepath<-"D:/varasto_iceagesimu"

hadmonths=1 hadoperaatio=1 ## 1: tex tjuly or precipjuly, 3 annual rain hadfilename="test.nc" hadbaseyear=0

1. hadvariaabeli="pr"

hadvariaabeli="tas"

unitti="Celcius" hadlonkkunimi="Near-Surface Air Temperature" hadlonvar1="longitude" hadlatvar1="latitude" hadmonth=fmonth hadyear=fage hadyears=numyears

haditem=15006 ## !!! will override in many cases!

## hadcm3b 60ka files path

hadbasepath<<-"D:/varasto_iceagesimu" hadbaseyear=-1 hadprocesspath<-"./data_processing/" lones1=0 latis1=0 inputdataname1=""

1. 1. non control vars

capture_hadcm3_temperature=0 capture_hadcm3_rainfall=0

1. 1. input data dirs
   2. etopo1 dir

predata="./predata/"

1. 1. Trace21ka TS dir

predata2="./predata2/"

dir.create("./data_processing") dir.create("./plotz")

1. outpath11<-"./data_processing/area.nc"
2. outpath12<-"./data_processing/area_neworog.nc"
   1. glac1D HDC file

inpath_tarasov1<-"./predata/TOPicemsk.GLACD26kN9894GE90227A6005GGrBgic.nc"

1. 1. peltier ice thickness 122 kyr

inpath_peltier1<-"./predata/IceT.I6F_C.131QB_VM5a_1deg.nc"

1. outpath22<-"./data_processing/oroin.nc"
   1. ! grid-registreted etopo1
2. inpath11<-"./predata/ETOPO1_Ice_g_gdal.grd"
   1. cell reg maybe inaccurate, but func w/dissever

inpath_etopo1<-"./predata/ETOPO1_Ice_c_gdal.grd" inpath_etopo2<- "./predata/etopo1_360.nc" data_processing="./data_processing/"

invarfname1<-"./data_processing/tempin.nc" invarfname2<-"./data_processing/precipin.nc"

inorofname2<-"./data_processing/area.nc" inorofname3<-"./data_processing/area_neworog.nc"

outorofname1<-"./data_processing/area_glacial.nc" outlandseamaskfname1<-"./data_processing/landsea_glacial.nc" outlandseamaskfname2<-"./data_processing/landsea_glacial.gif"

outvarfname1<-"./data_result/result.nc" outvarfname1b<-"./data_result/result_land.nc" outvarfname2<-"./data_result/result2.nc"

outvarfname_accurate1<-"./data_result/accurate.nc"

inorofname360_1<-"./data_processing/oroin.nc" invarfname360_1<-"./data_processing/tempin.nc" inicefname360_1<-"./data_processing/maskin.nc"

inorofname180_1<-"./data_processing/oroin_180.nc" invarfname180_1<-"./data_processing/tempin_180.nc" inicefname180_1<-"./data_processing/maskin_180.nc"

invarfname180_2<-"./data_processing/precipin_180.nc"

smallrainame1<-"./data_processing/smallrain.nc"

1. outname1<-"./data_processing/chelsa_current_annual_rain.nc"
2. outname2<-"./data_processing/chelsa_lgm_ccsm4_annual_rain.nc"

outname3<-"./data_processing/etopo1.nc" outname4<-"./data_processing/lons.nc" outname5<-"./data_processing/lats.nc" outname6<-"./data_processing/distance.nc" outname7<-"./data_processing/slope.nc" outname8<-"./data_processing/aspect.nc" outname9<-"./data_processing/hillshade.nc" outname10<-"./data_processing/tpi.nc" outname11<-"./data_processing/westness.nc" outname12<-"./data_processing/blurelev.nc" outname13<-"./data_processing/windir.nc" outname14<-"./data_processing/etopo_big.nc"

plottauzdir1="./plotz/" plotfname0="result0.svg" plotfname1="./plotz/result1.svg"

resultdir1<-"./data_result/"

dir.create(data_processing) dir.create(plottauzdir1) dir.create(resultdir1)

1. 1. csv related variables, not yet used

smallrst="./data/small/" desticsv="./data/" destirst="./data/rst/" smalltemperaturepath<-"./data/small/small_temperature.nc" elevpath<-"./data/rst/elevation.nc" smallelevpath<-"./data/small/small_elevation.nc" distancepath<-"./data/rst/distance.nc" smalldistancepath<-"./data/small/small_distance.nc" latpath<-"./data/rst/latitude.nc" smalllatpath<-"./data/small/small_latitude.nc" lonpath<-"./data/rst/longitude.nc" smalllonpath<-"./data/small/small_longitude.nc"

1. 1. CODE BEGINS

exte1<- extent(lon1,lon2,lat1,lat2)

1. 1. global doordinates: pacific or europe centered
   2. rasnum 1: -180 - 180
   3. rasnum 2: 0 - 360

rasnum1=1

if(lon1>180)

  {

rasnum1=2;

  }
   
if(lon2>180)
   {

rasnum1=2; }

1. 1. 1. 1. 1. 1. 1. 1. 1. 2. subroutines

sealevel_from_age<-function(inage) {

1. 1. agetable1<-c(60000,24000,22000,20000,18000,16000, 15000,14000,13000,12000,11000, 10000,8000,7000,6000 ,0)
   2. leveltable1<-c(-65,-125,-130,-125,-120,-115, -110,-80,-75,-65,-60, -45,-15,-10,-5 ,0 )

agetable1<-c(110000, 65000,60000,55000,50000,45000, 40000,37500, 35000, 30000, 24000,22000,20000,18000,16000, 15000,14000,13000,12000,11000, 10000,8000,7000,6000 ,0) leveltable1<-c(-40, -85,-75, -60,-70,-80,-75, -65, -80, -75, -125,-130,-125,-120,-115, -110,-80,-75,-65,-60, -45,-15,-10,-5 ,0 )

1. plot(agetable1, leveltable1, main = "Sea level curve")

#points(approx(agetable1, leveltable1), col = 2, pch = "*")

# 	selev1=points(approx(agetable1, leveltable1), col = 2, pch = "*")

seali<<-approx(agetable1, leveltable1, xout=inage)

    sealevel<-as.numeric(seali)

   print("Calculated sea level ")
   print(seali)
   print(sealevel)

}

default_settings_on<-function() { print("Default settings on.") arg_reading <<-1 orography_preprocess<<-1 capture_temperature<<-1 capture_elevation<<-1 normal_ds<<-1 global_ds<<-0 accurate_ds<<-0 plot_result<<-1 }

delete_directories<-function() { print("Delete dirs") getwd() unlink(data_processing, recursive = TRUE) unlink(plottauzdir1, recursive = TRUE) unlink(resultdir1, recursive = TRUE)

}

read_argus<-function() {

if (arg_reading >0) { args = commandArgs(trailingOnly=TRUE)

if (length(args)==0) { ## stop("Not args. "=FALSE) print("usage like: rscript --vanilla tracs3.r 10000") # stop("Abort.") }

if (length(args)>0) {

print("Reading of argus done.")

eka<-args[1] toka<-args[2] fage<<-as.numeric(eka) print(fage)

delete_directories() default_settings_on() #fmonth=as.numeric(toka) } }

} ## argus

1. 1. 1. MAIN PROGE

writeout<-function(oras, outn, varnamex, varunitx, longnamex) {

crs(oras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(oras, filename=outn, overwrite=TRUE, format="CDF", varname=varnamex, varunit=varunitx, zname="Zname", longname=longnamex, xname="lon", yname="lat")

}

preprocess_orography<-function() {

print("Processing orography, wait ...") rin1<-raster(inpath_etopo1)

rin1

x1 <- crop(rin1, extent(-180, 0, -90, 90)) x2 <- crop(rin1, extent(0, 180, -90, 90)) extent(x1) <- c(180, 360, -90, 90) rin2 <- merge(x1, x2)

rin2

    rasnum1=1
    
    if(lon1>180)
    {

rasnum1=2; }

      if(lon2>180)
    {

rasnum1=2; }

   if(rasnum1==1)

{ rout1<- crop(rin1, exte1) }

else { rout1<- crop(rin2, exte1) }

rout1

1. plot(rout1)

   heightdelta=sealevel*-1
   

rout2=rout1+heightdelta

rout2[rout2 < 0] <- 0

plot(rout2, col=viridis(100))

1. inorofname2<-"./data_processing/area.nc"
2. inorofname3<-"./data_processing/area_neworog.nc"

crs(rout2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rout2, inorofname3, overwrite=TRUE, format="CDF", varname="Elev", varunit="meters",

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

crs(rin2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rin2, inpath_etopo2, overwrite=TRUE, format="CDF", varname="Elev", varunit="meters",

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

writeRaster(rout1, inorofname2, overwrite=TRUE, format="CDF", varname="Elev", varunit="meters",

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

1. Pazka

# rout2 #plot(rout2)

}

generate_hadfilename<-function(hadbaspath1, yrr1, varr1) { #hadfilename="D:/varasto_iceagesimu/bias_regrid_tas_55_57.5kyr.nc" hadfilenamex1=hadbaspath1

hadfilenamex1<-paste0(hadfilenamex1,"/bias_regrid_") hadfilenamex1<-paste0(hadfilenamex1,varr1) hadfilenamex1<-paste0(hadfilenamex1,"_")

yrr1=yrr1-(2500/2)+1 a=yrr1/100 b=a c=round(b/25,0) d=(c*25)/10 e=d+2.5 hadbaseyear<<-d*1000

#print (b) #print (c) #print (d) #print (e)

hadfilenamex1<-paste0(hadfilenamex1,toString(d)) hadfilenamex1<-paste0(hadfilenamex1,"_") hadfilenamex1<-paste0(hadfilenamex1,toString(e)) hadfilenamex1<-paste0(hadfilenamex1,"kyr.nc") return(hadfilenamex1) }

1. 1. loadvar 3d

load_trace_var3d_3<-function(ivars1, lons1, lats1,variaabeli, unitti, age1, fage,numyears, fmonth1)

{

1. ivars1 <- ncvar_get(input1, variaabeli)
2. lats1 <- ncvar_get(input1, "lat")
3. lons1 <- ncvar_get(input1, "lon")

print("Input variable processing ...")

  	varlocation1<<-((fage-age1)*12*-1)+fmonth1
  

if (numyears==1) { vaari1 <-ivars1[,,varlocation1] }

if (numyears>1) { ## vaari10<-ivars1[,,varlocation1]

vaari10=0 for (nn in 0:(numyears-1)) { vaari10 <-vaari10+ivars1[,,varlocation1+nn*12] }

vaari1=vaari10/numyears }

   if (variaabeli=="TS")
   {

vaari1<-vaari1-273.15

    }
    

vaari1<- apply(t(vaari1),2,rev)

rrvar1<-raster (vaari1)

rrvar1@extent<-extent(0, 360, -90, 90)

  print("Write ..")

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

writeRaster(rrvar1, invarfname360_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

   rrvar1_180<-rotate(rrvar1)

   rrvar1_180@extent<-extent(-180, 180, -90, 90)

plot(rrvar1_180, col=rainbow(120))

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

writeRaster(rrvar1_180, invarfname180_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

#stop("Joo")

}

tracelocation3 <- function(dirra1, variaabeli, fage, numyears, fmonth) {

   aage2=fage-numyears
   aage1=fage
   

lista0<-list.files(path=dirra1,pattern="trace*", ignore.case=TRUE) pitu<-length(lista0)

1. print(fage1,)

#print("Kup")

   #print (pitu)
   
#     lista0

fagek1=as.integer(aage1) fagek2=as.integer(aage2)

if(fagek2>fagek1) { huba=fagek1 fagek1=fagek2 fagek2=huba

}

mara1=-999 mara2=-999

for (n in 1:(pitu)) { row1<-lista0[n] jono1<-row1 #print (jono1)

     	s1<-gsub('\\.', '_', jono1[1])

s2<-gsub('\\-', '_', s1) s3<-gsub('BP', , s2)

q1 <-str_split(s3, "_")

      	sr1<-q11

age10<-as.integer(sr1[3])

		age20<-as.integer(sr1[4])
		
		age1<-age10[1]
		age2<-age20[1]
		                  
         if (fagek2==400) 

{ fagek2<-0 }

if(fagek1<age1)

    		{
			if(fagek1>age2)	

{ mara1<-n } }

if(fagek1==age1)

    		{

mara1<-n # break }

if(fagek2<age1) { if(fagek2>age2) { mara2<-n break } }

if(fagek2==age1)

    		{

mara2<-n }

1. print (paste0(str(age1)," ",str(age2)))

}

print("Loop done")

print(mara1) print(mara2)

1. Haista

ivars0=0 ivars1=0 mara=0 natta=0

agari1=0

1. 1. jn warning debug

for (mara in mara1:mara2)
  {

row1<-lista0[mara] print(row1)

jono1<-row1 s1<-gsub('\\.', '_', jono1[1]) s2<-gsub('\\-', '_', s1) s3<-gsub('BP', , s2) q1 <-str_split(s3, "_") sr1<-q11

age1<-as.numeric(sr1[3]) age2<-as.numeric(sr1[4])

# print (mara) filename1<-paste0(predata2, row1) print(filename1)

    	putin1 <- nc_open(filename1)
    	
    	

print(class(putin1))

    #     stop("k")
   #   ivars0<-0

ivars0 <- ncvar_get(putin1, variaabeli) lones1<- ncvar_get(putin1, lonvar1) latis1<- ncvar_get(putin1, latvar1)

if(natta>0) { ivars1<-abind(ivars1,ivars0,along=3) #ivars1=c(ivars1,ivars0) } else {

agari1<-as.numeric(age1) ivars1<-ivars0 }

       natta=natta+1
 
       nc_close(putin1)
  }	

elevlocation1<<-round( (fage-26000)/100)*-1

  class(ivars1)
  dim (ivars1)
  
  #	varlocation1<<-((fage-age1)*12*-1)+fmonth

1. 1. orig
2. load_trace_var3d_3(ivars1, lons1, lats1,variaabeli, unitti,agari1, fage, numyears, fmonth )

load_trace_var3d_3(ivars1, lons1, lats1,variaabeli, unitti,agari1, fage ,numyears, fmonth)

   }

load_glac_elev<-function() {

number2<-elevlocation1

input2 <- nc_open(inpath_tarasov1)

class(input2)

#input2

print("INPUT GLAC1D ELEV")

elev <- ncvar_get(input2, "HDC") lats2 <- ncvar_get(input2, "YLATGLOBP5") lons2 <- ncvar_get(input2, "XLONGLOB1")

print("INPUT GLAC ELEV DONE")

   print(number2)

elev0<-elev[,,number2]

str(elev0)

elev1<-apply(t(elev0),2,rev) elev1[elev1 <0.0] <- 0.0

rr2<-raster(elev1)

rr2@extent<-extent(0, 360, -90, 90)

2. plot(rr2, col=rainbow(64))

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

writeRaster(rr2, inorofname360_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

   rr2_180<-rotate(rr2)
   
   rr2_180@extent<-extent(-180, 180, -90, 90)

plot(rr2_180, col=rainbow(120))

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

writeRaster(rr2_180, inorofname180_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

}

1. 1. load_var_nc

load_glac_icemask<-function() {

number2<-elevlocation1

1. inpath22

input2 <- nc_open(inpath_tarasov1)

1. class(input2)

#input2

print("INPUT GLAC1D ICEMASK")

elev <- ncvar_get(input2, "ICEM") lats2 <- ncvar_get(input2, "YLATGLOBP5") lons2 <- ncvar_get(input2, "XLONGLOB1")

print("INPUT GLAC ICEMASK DONE")

elev0<-elev[,,number2] elev1<-apply(t(elev0),2,rev) elev1[elev1 <0.0] <- 0.0

rr2<-raster(elev1)

rr2@extent<-extent(0, 360, -90, 90)

2. plot(rr2, col=rainbow(64))

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

writeRaster(rr2, inicefname360_1, overwrite=TRUE, format="CDF", varname="Icem", varunit="percent", longname="Icemask", xname="lon", yname="lat")

   rr2_180<-rotate(rr2)
   
   rr2_180@extent<-extent(-180, 180, -90, 90)

plot(rr2_180, col=rainbow(120))

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

writeRaster(rr2_180, inicefname180_1, overwrite=TRUE, format="CDF", varname="Icem", varunit="percent", longname="Icemask", xname="lon", yname="lat")

}

1. 1. load_var_nc

load_peltier_elev<-function() { print("Input Peltier 6G 122 elevation data.")

input2 <- nc_open(inpath_peltier1)

ele <- ncvar_get(input2, "stgit") lats2 <- ncvar_get(input2, "Lat") lons2 <- ncvar_get(input2, "Lon") time2 <- ncvar_get(input2, "Time")

   possi2<- -999
   leenu2<-length(time2)
   hage<-(fage/1000)
     

for (n2 in 1:(leenu2) ) { taa2<-time2[n2]

if(hage==taa2) { possi2=n2 break }

if(hage>taa2) { possi2=n2-1 break }

}

elev0<-ele[,,possi2]

elev1<-apply(t(elev0),2,rev) elev1[elev1 <0.0] <- 0.0

rr2<-raster(elev1)

rr2@extent<-extent(0, 360, -90, 90)

plot(rr2, col=rainbow(64))

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

writeRaster(rr2, inorofname360_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

   rr2_180<-rotate(rr2)
   
   rr2_180@extent<-extent(-180, 180, -90, 90)

plot(rr2_180, col=rainbow(120))

crs(rr2_180) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 1 writeRaster(rr2_180, inorofname180_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

## peltier ice6g elevation

}

nok_process_rasters<-function() { # WARNING!!!! vprocedure is NOT OK

## NOTE DISTANCE CALCULATION IS OFF SAKE FOR SPEED ## if you use distance, you must it below in script

re_in <- raster("./data_processing/oroin.nc") rt_in<-raster("./data_processing/tempin.nc")

re_in[re_in <0.0] <- 0.0

plot(re_in) plot(rt_in)

#extent(re_in)<-(0,360,-90,90)

big_sabluna<- raster(nrow=180, ncol=360)

#stop("TEST BRK")

extent(big_sabluna)<-c(0,360,-90,90) res(big_sabluna)<-c(1.0, 1.0)

#small_sabluna<-raster(nrow=48, ncol=96) #extent(small_sabluna)<-c(0, 360.000, -90, 90 ) #res(small_sabluna)<-c(3.75, 3,78947) #small_sabluna<-raster(nrow=50, ncol=100) #small_sabluna<-raster(nrow=5, ncol=10) small_sabluna<-raster(nrow=45, ncol=90)

extent(small_sabluna)<-c(0, 360.000, -90, 90 )

1. extent(small_sabluna)<-c(-1.875, 358.125, -89.01354, 89.01354 )
2. res(small_sabluna)<-c(3.75, 3.708898)

re_big <- resample(re_in,big_sabluna, method='bilinear') extent(re_big)<-c(0,360,-90,90) rt_big <- resample(rt_in,big_sabluna, method='bilinear')

1. extent(re_big)<-c(0,360,-90,90)

re_small<-resample(re_in,small_sabluna, method='bilinear')

extent(re_small)<-c(0, 360, -90, 90 )

1. res(re_small)<-c(3.75, 3.75)

rt_small <- resample(rt_in,small_sabluna, method='bilinear')

extent(rt_small)<-c(0,360,-90,90)

1. res(rt_small)<-c(3.75, 3.75)

rlon_big <- rlat_big<- re_big xy <- coordinates(re_big) rlon_big[] <- xy[, 1]

1. JN warning

rlat_big[] <- xy[, 2]

rlat_small<-resample(rlat_big,small_sabluna, method='bilinear') rlon_small<-resample(rlon_big,small_sabluna, method='bilinear')

1. plot(rlat_small)

1. re_big<-rotate(re_big)

1. extent (re_big)

1. plot(re_big)

1. plot (rb_out)
2. extent(rb_out)

1. extent(rt_in)
2. rt_in

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

writeRaster(rt_small, smalltemperaturepath, overwrite=TRUE, format="CDF", varname="tempa", varunit="Units",

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

1. plot (rt_in)

plot (rt_small)

rt_in

rt_small

1. re_big<-rotate(re_big)
2. re_small<-rotate(re_small)
3. rt_big<-rotate(rt_big)
4. rt_small<-rotate(rt_small)
5. rlon_big<-rotate(rlon_big)
6. rlon_small<-rotate(rlon_small)
7. rlat_big<-rotate(rlat_big)
8. rlat_small<-rotate(rlat_small)

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

writeRaster(re_big, elevpath, overwrite=TRUE, format="CDF", varname="elev", varunit="Units",

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

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

writeRaster(re_small,smallelevpath, overwrite=TRUE, format="CDF", varname="elev", varunit="Units",

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

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

writeRaster(rlat_big, latpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

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

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

writeRaster(rlat_small,smalllatpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

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

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

writeRaster(rlon_big, lonpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

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

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

writeRaster(rlon_small, smalllonpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

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

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

writeRaster(rt_big, temperaturepath, overwrite=TRUE, format="CDF", varname="tempa", varunit="Units",

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

r_1=re_big r_1[r_1 > 0] <- NA r_1[r_1 < 1] <- 1

if(kalk_distance==1) { rdistance_big <- distance(r_1)

# meters to degrees

rdistance_big<-rdistance_big/111120

rdistance_small<-resample(rdistance_big,small_sabluna, method='bilinear')

rdistance_big<-rotate(rdistance_big) rdistance_small<-rotate(rdistance_small)

#plot(rdistance_big)

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

writeRaster(rdistance_small, smalldistancepath, overwrite=TRUE, format="CDF", varname="distance", varunit="Units", longname="Distance", xname="lon", yname="lat")

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

writeRaster(rdistance_big, distancepath, overwrite=TRUE, format="CDF", varname="distance", varunit="Units", longname="Distance", xname="lon", yname="lat") }

} ## process_rasters

nok_calculate_tables<-function() {

# WARNING procedure is NOT OK

print ("Tables.")

climdir="./data/" smalldir="./data/small/"

elevras_small<-raster("./data/small/small_elevation.nc") distras_small<-raster("./data/small/small_distance.nc") tempras_small<-raster("./data/small/small_temperature.nc")

#elevras_small #distras_small #tempras_small

temp_df<-as.data.frame(tempras_small, xy = TRUE) # jn warn

temp.df[is.na(temp.df)] <- 0.0 elev_df<-as.data.frame(elevras_small, xy = TRUE) dist_df<-as.data.frame(distras_small, xy = TRUE)

#temp_df

plot(tempras_small, col=rainbow(100)) #plot(elevras_small) #plot(distras_small)

small_stack <-stack(elevras_small, tempras_small, distras_small)

#plot(small_stack)

small_df<-data.frame( rasterToPoints( small_stack ) )

#plot(small_df)

## JN warning !!!!

small.df[is.na(small.df)] <- 0.0

#small_df

#colnames(small_df) <- c(library(reshape2) colnames(small_df) <- c("Lon","Lat","Elev", "AvgTemp","distCoast" )

#temp_df

#plot(tempras_small)

colnames(small_df) <- c("Lon","Lat","Elev", "AvgTemp","distCoast" )

small_df["StationID"] <- seq(1,nrow(small_df))+10

#small_df2<- small_df[c("Lat","Lon","Elev", "AvgTemp","distCoast" )]

#StationName,StationID,Lat,Lon,Elev,AvgTemp,distCoast

#colnames(small_df2 )

#head(small_df2,5)

small_df2 <- small_df[c("StationID","Lon","Lat","Elev", "AvgTemp","distCoast" )]

#small_df2

#len(small_df2.AvgTemp) #len(small_df2.Lon)

colnames(small_df2 )

head(small_df2,5)

sample_df=sample_n(small_df2, 600)

## JN warning !!!!

sample_df[is.na(sample_df)] <- 0.0

write.csv(small_df2,"./data/full.csv") write.csv(sample_df,"./data/sample.csv")

##rekt<-dplyr::filter(mag <5.5 & mag > 4.5)

#dplyr::filter(small_df2), lat>

d1<-filter(small_df, Lon >-15.0) d2<-filter(d1, Lon <40.0) d3<-filter(d2, Lat <70.0)

europedf<-filter(d3, Lat >30.0)

write.csv(europedf,"./data/europe.csv")

#plot(europedf)

}

1. 1. Downscaling utilities

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

   print ("Downscaler()")			

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)

min_iter <- 5 # Minimum number of iterations max_iter <- 10 # Maximum number of iterations p_train <- 0.1 # Subsampling of the initial data

	 if(method>9999)
	 {

method=2 }

## dissever run

   if(method==2)

{ oma_juttu <- dissever(coarse = coarse_raster, fine = raster_stack, method = "lm", p = p_train, min_iter = min_iter,max_iter = max_iter, verbose=1) orotemp<-oma_juttu$map }

## spatialeco downscale if(method==1) { oma_juttu2 <- raster.downscale(p1, coarse_raster) orotemp<-oma_juttu2$downscale }

1. 1. delta regression 1,1

if(method==0) {

orotemp<-orodelta

   	}

1. 1. delta regression by lapse rate

if(method==3) { orotemp<-orodelta*0.0065*-1

   	}

#biassi=4

#tempiso<-baset1+oma_juttu$map+biassi

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

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

1. plot(outtemp, col=rev(rainbow(256)) )

outtempr<-rotate(outtemp)

plot(outtempr)

     return(outtemp)

}

downscale_orography<-function(method) {

  print("Orography ..")

1. 1. downscale orography

   oroko1<<-inorofname180_1;
   
   rasnumi1<<-rasnum1    
   
   "rasnum1"
   print(rasnumi1)
    
   
    if(rasnumi1==2)
    {

print ("ORO 0-360") oroko1<<-inorofname360_1; }

coarse_raster1<-raster(oroko1) fine_raster1<-raster(inorofname3)

if(capture_elevation==2) { print("Peltier test code") vali_raster1<-resample( fine_raster1, coarse_raster1) coarse_raster1<-coarse_raster1+vali_raster1 }

1. coarseoro<- resample(p1, coarse_raster)

  kover_raster1<-resample(coarse_raster1, fine_raster1)
  kover_koeff1<-kover_raster1
  
  kover_koeff1[kover_koeff1>0]<- 1

   plot(fine_raster1)
   plot(coarse_raster1)

   names(fine_raster1)<-"Oroa"

## because europe lies between -15 ... 40, must rotate

1. 1. if(rasnumi1==1)

##    {

1. 1. "rotate oro"
   2. coarse_raster1<-rotate(coarse_raster1)
   3. }

    coarse_raster1[coarse_raster1 < 0] <- 0

plot(fine_raster1, col=viridis(100)) plot(coarse_raster1, col=viridis(100))

outras1<-downscaler(coarse_raster1, fine_raster1,method)

   outras1[outras1<1]<- 0

#"Oro ds debug 1"

1. names(outras1)<-"Oro"

1. sabluna1<-raster(inorofname3)

 #"sabluna1"
#  plot(sabluna1)
  
#  sabluna1[sabluna1<0]<- 0

1. sabluna1[sabluna1>0]<- 1

1. outras1<-outras1*sabluna1

#  outras1<-outras1*kover_koeff1
  
#  plot(outras1)
  
 

1. 1. elevenaame1<-"./data_processing/area_neworog.nc"
2. elevenaame2<-"./data_processing/area_glacial.nc"

plot(outras1, col=rev(plasma(256)))

writeRaster(outras1, outorofname1, overwrite=TRUE, format="CDF", varname="Oro", varunit="meters", longname="Orogr", xname="lon", yname="lat")

   outras2<-outras1
   
 
   outras2[outras2[]<=0.0] <- 0 
   outras2[outras2[]>0.0] <- 1.0 
 	plot(outras2, col=rev(plasma(256)))
 	
 	writeRaster(outras2, outlandseamaskfname1, overwrite=TRUE, format="CDF", varname="Land", varunit="unit", 

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

1. writeRaster(outras2, outlandseamaskfname2, overwrite=TRUE, format="CDF", varname="Land", varunit="unit",
2. longname="Land", xname="lon", yname="lat")

}

downscale_temperature<-function(method) {

1. 1. downscale temperature

   print("DS of variable...")
   
    rasnumi1<<-rasnum1   

varoko1<<-invarfname360_1;

if(rasnumi1==1) {

1. coarse_raster2<-rotate(coarse_raster2)

        varoko1<<-invarfname180_1;

}

coarse_raster2<-raster(varoko1)

fine_raster2<-raster(outorofname1)

print("TS rsaters")

#coarse_raster2

#fine_raster2

1. if(klopaali1==1)
2. {
3. fine_raster2=rotate(fine_raster2)
4. }

plot(fine_raster2)

   plot(coarse_raster2)

#fine_raster2<-outras1

1. fine_raster2<-raster(inorofname2)
   1. if(rasnumi1==1)

##    {

1. 1. coarse_raster2<-rotate(coarse_raster2)
   2. }

plot(fine_raster2, col=viridis(100)) plot(coarse_raster2, col=viridis(100) )

outras2<-downscaler(coarse_raster2, fine_raster2,method)

1. names(outras1)<-"TS"

plot(outras2, col=rev(rainbow(256)))

1. writeRaster(outras2, filename='./outdata/result.nc', format="CDF", overwrite=TRUE)

1. out_text_1="TS"
2. out_text_2="oC"
3. out_text_3="Temperature"

writeout(outras2,outvarfname1,out_text_1, out_text_2, out_text_3)

mask1<-raster(outlandseamaskfname1) outras3<-outras2 outras3<-outras3*mask1 writeout(outras3,outvarfname2,out_text_1, out_text_2, out_text_3)

#writeRaster(outras3,outvarfname2 , overwrite=TRUE, format="GIF", varname="TS", varunit="Celcius", #longname="Temperature", xname="lon", yname="lat")

}

normal_ts_downscale<-function()

     {

## normal ## methods 0 delta, 1 spatialeco 2 dissever 3 tempdelta

print("Normal etopo1 ds. Orography ...")

  if(capture_elevation==2)
  {	

# force method to 0 because of peltier topography downscale_orography(0)

  }
  else
  {	

downscale_orography(method_ds_oro)

  }

print("Temperature ...") downscale_temperature(method_ds_var)

1. 1. spatialeco
   2. downscale_orography(1)
   3. downscale_temperature(1)

1. 1. delta

1. 1. downscale_orography(0)
2. downscale_temperature(3)

}

world_ts_downscale<-function()

     {

rasnum1<<- 2 klopaali1<<-1

inorofname1<<-"./data_processing/oroin.nc" invarfname1<<-"./data_processing/tempin.nc"

1. 1. inorofname1<-"./data/rst/elevation.nc"
   2. invarfname1<-"./data/small/small_temperature.nc"

## normal print("WORLD DS") # downscale_orography(2) ## WORLD DS #outorofname1<<-"./data/rst/elevation.nc" outorofname1<<-"./data_processing/oroin.nc" # outorofname2<<-"./data_processing/oroin.nc"

downscale_temperature(3) }

custom_ts_downscale_1<-function(orofilee1, tempfilee1)

     {

print("Custom ds type 1.") print(orofilee1) print(tempfilee1)

rasnum1<<- 1

    coarse_raster=raster(tempfilee1)
    fine_raster=raster(orofilee1)
   	  
    method=3

outras<-downscaler(coarse_raster, fine_raster,method)

plot(outras, col=rev(rainbow(256)))

writeRaster(outras, "./data_result/accurate.nc", overwrite=TRUE, format="CDF", varname="TS", varunit="Celcius", longname="Temperature", xname="lon", yname="lat")

}

plottaus1<-function(zvarnaame1, outfilenaame1) {

1. elevenaame0<-"./data_processing/area_neworog.nc"

elevenaame1<-"./data_processing/area_glacial.nc"

   while (dev.cur()>1) dev.off()

grayscale_colors <- gray.colors(100,start = 0.0,end = 1.0,gamma = 2.2,alpha = NULL)

relev1<-raster(elevenaame1)

razteri1<-raster(zvarnaame1)

1. plot(relev1)

aage<<-fage

ram1<-relev1

maintitle1<-paste0("Temperature of July, ",toString(aage), " BP") subtitle1<-"Trace21ka CCSM3 downscaled" xlab1<-"Lon" ylab1<-"Lat"

print (maintitle1)

limx1<- -15 limx2<-45

   limy1<-35
   limy2<-65	
      
   limitx1<-c(limx1, limx2)
   limity1<-c(limy1, limy2)
   
   templevs1<-c(-15,-10,-5,0,5,10,15,20,25,30,35)
   ltemp1<-length(templevs1)

tempmin1<-templevs1[0]

   tempmax1<-templevs1[ltemp1-1]

print(tempmin1) print(tempmax1)

x1 <- rasterToContour(razteri1, levels=templevs1)

1. 1. lev1 <- seq(tempmin1,tempmax1, by=5)

 lev1=seq(-15,35,5)

class(x1) ek1<-rasterToContour(ram1)

svg(filename="out.svg", width=12, height=10, pointsize=20)

plot( razteri1, main=maintitle1, sub=subtitle1, xlab=xlab1, ylab=ylab1, xlim=limitx1, ylim=limity1, breaks=templevs1, col=rev(rainbow(ltemp1)) )

#  plot(x1,xlim=c( -15, 50 ), ylim=c( 30, 60 ), alpha=0.5,add=TRUE)

plot(ram1, xlim=limitx1, ylim=limity1 , zlim=c(0,1), breaks=c(-1,0,1), col=viridis(2) , legend=FALSE, add=TRUE)

1. plot(ek1, add=TRUE, legend=FALSE)

dev.off() }

fill.na <- function(x) { width1<-67

 center = 0.5 + (width1*width1/2) 
 if( is.na(x)[center] ) {
   return( round(mean(x, na.rm=TRUE),0) )
 } else {
   return( round(x[center],0) )
 }

}

get_europe_gtopo30<-function() { # france #ext1 <- extent(-4,8, 39, 51) #dordogne #ext1 <- extent(-2,4,42 , 48)

#europe ext1 <- extent(-15,40,30,70)

outname1<-"./predata/europe30.nc"

r1<-raster("./predata/gt30w020n90.tif")# r2<-raster("./predata/gt30e020n90.tif") # r4<-raster("./predata/gt30w020n40.tif") # r3<-raster("./predata/gt30e020n40.tif")

#plot(r1)

mosaik1 <- merge(r1,r2,r3,r4)

kropped1<-crop(mosaik1, ext1)

#plot(kropped1)

crs(kropped1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(kropped1, filename=outname1, format="CDF", overwrite=TRUE) }

load_hadcm3_slice <- function(hadfilename, hadvariaabeli, haditem) { print("Trying to load hadcm3 filee ")

## varlocation1<<-((fage-age1)*12*-1)+fmonth1

varlocation1=haditem

print(hadfilename)

hadcount=1

    	putin1 <- nc_open(hadfilename)
    	  		

print(class(putin1))

lones1<- ncvar_get(putin1, hadlonvar1) latis1<- ncvar_get(putin1, hadlatvar1) t <- ncvar_get(putin1, "time")

lenlones1<-length(lones1) lenlatis1<-length(latis1)

vaari0<- ncvar_get(putin1,hadvariaabeli, start=c(1,1,haditem), count=c(lenlones1,lenlatis1,1) )

nc_close(putin1)

print(dim(vaari0))

padding1 = matrix(0,lenlones1,180)

vaari1<-cbind(padding1,vaari0) #vaari1<-vaari0

print(dim(vaari1))

  vaari1<- apply(t(vaari1),2,rev)
  

rrvar00<-raster (vaari1)

rrvar00@extent<-extent(0, 360, -90, 90)

width1 = 67 rrvar1<- focal(rrvar00, w = matrix(1,width1,width1), fun = fill.na,

           pad = TRUE, na.rm = FALSE)
           
           

summary(getValues(rrvar1))

plot(rrvar1, col=rainbow(120))

  print("invarfname:")
  print(invarfname360_1)
  
  print("Write hadcm3 file ..")

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

writeRaster(rrvar1, invarfname360_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

   rrvar1_180<-rotate(rrvar1)

   rrvar1_180@extent<-extent(-180, 180, -90, 90)

plot(rrvar1_180, col=rainbow(120))

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

writeRaster(rrvar1_180, invarfname180_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

  # print("Test break")
   
  # stop()

   }

load_hadcm3_year_month <- function(hadfilename, hadvariaabeli, hadyear, hadmonth, hadbeginyear) {

hadyeardelta=hadyear-hadbeginyear haddelta=hadyeardelta*12+hadmonth haditem=30001-haddelta

load_hadcm3_slice(hadfilename, hadvariaabeli, haditem) }

load_average_hadcm3_years_months <- function(hadfilename, hadvariaabeli,hadbeginyear, hadyear, hadmonth, hadyears, hadmonths, operaatio, metoden1) {

## operaatio==0: sum ## operaatio==1: averege of data matrices

print ("Loading many HADCM3 slices.") print(hadfilename)

	putin1 <- nc_open(hadfilename)
    	  		          

lones1<- ncvar_get(putin1, hadlonvar1) latis1<- ncvar_get(putin1, hadlatvar1) t <- ncvar_get(putin1, "time") lenlones1<-length(lones1) lenlatis1<-length(latis1)

plastic1<-matrix(0,lenlones1,lenlatis1)

kohtia=0

hadyearend=hadyear+hadyears hadmonthend=hadmonth+hadmonths

#print ("Hadyear") #print (hadyear) #print (hadyearend) #print (hadmonth) #print (hadmonthend)

hadyeardelta1=hadyear-hadbeginyear haddelta1=hadyeardelta1*12 haditem1=30001-haddelta1

hadpususize1=hadyears*12

#print ("Haddeltas ") #print( hadyeardelta1) #print(haddelta1) #print(haditem1)

print("Get had data ...")

hadpusu1<-ncvar_get(putin1,hadvariaabeli, start=c(1,1,haditem1), count=c(lenlones1,lenlatis1,hadpususize1) )

   nc_close(putin1)	 

#print (hadpusu1[7])

hadex1=0

   ammax=hadyears-1
   kammax=hadmonth+hadmonths-1

#print("Extract haddata:")

hadpusupitu1=length(hadpusu1)

#print (hadpusupitu1) #print (hadpususize1)

#print (ammax) #print (kammax) #print (hadmonth)

#print ("luup")

vuozia1=0

for (m in (1:ammax) ) { #print("R") for (n in (hadmonth:kammax) ) { #print(".")

#print(n)

hadex1=m*12+n

# print (hadex1) vaari0<-hadpusu1[,,hadex1] plastic1<-plastic1+vaari0

kohtia=kohtia+1

} vuozia1=vuozia1+1 }

print("Kohtia:") print(kohtia) print("Vuozia1:") print (vuozia1) print("Operaatio:") print(operaatio)

if(operaatio==0) { ## operaatio 0, sum ## one year precip print("sum") vaari3<-plastic1 }

if(operaatio==1) { print("one month, many years") ## one month, many years vaari3<-(plastic1/kohtia) # JN WARNING koe }

if(operaatio==2) { print("annualaverage, one year") ## annual #print("Operaatio 3: temperature annual") vaari3<-plastic1/12 }

if(operaatio==3) { print("annual sum, many years") #print("Operaatio 3: rainfall annual") vaari3<-plastic1/vuozia1 }

#  plot(vaari3)

vaari4<-vaari3 vaari5<-vaari3 vaari6<-vaari3 vaari7<-vaari3

   dima1<-dim(vaari3)
   w1=dima1[1]
   h1=dima1[2]
   print (w1)
   print (h1)
  
  if(metoden1==2)
  { 
   

for (ix in (1:w1) )

{ beg1=1 viksu1=-1 for (iy in (1:h1) )

{ piksu1=vaari3[ix,iy]

if(is.na(piksu1)) {

if(beg1==0) { vaari4[ix,iy]=viksu1 vaari3[ix,iy]=1 } else { vaari4[ix,iy]=piksu1

}

} else { beg1=0 viksu1=piksu1 beg1=0 }

}

beg1=1 viksu1=-1 for (iy in seq(h1, 1, by=-1 ) ) {

piksu1=vaari3[ix,iy]

if(is.na(piksu1)) {

vaari4[ix,iy]=viksu1

} else ## ei na { beg1=0 viksu1=piksu1 beg1=0 }

} }

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. sekond

  				for (iy in (1:h1) )
   

{ beg1=1 viksu1=-1 for (ix in (1:w1) )

{ piksu1=vaari5[ix,iy]

if(is.na(piksu1)) {

if(beg1==0) { vaari6[ix,iy]=viksu1 vaari5[ix,iy]=1 } else { vaari6[ix,iy]=piksu1

}

} else { beg1=0 viksu1=piksu1 beg1=0 }

}

beg1=1 viksu1=-1 for (ix in seq(w1, 1, by=-1 ) ) {

piksu1=vaari5[ix,iy]

if(is.na(piksu1)) {

vaari6[ix,iy]=viksu1

} else ## ei na { beg1=0 viksu1=piksu1 beg1=0 }

} }

vaari7<-(vaari4+vaari6)/2 }

1. plot(vaari3)

   #print(dim(vaari3))

   padding1 = matrix(0,lenlones1,180)

   vaari1<-cbind(padding1,vaari7)  
   #vaari1<-vaari0

   #print(dim(vaari1))

   vaari1<- apply(t(vaari1),2,rev)
  

rrvar00<-raster (vaari1)

rrvar00@extent<-extent(0, 360, -90, 90)

1. 1. fill nas with nearest non-na

if(metoden1==1) { ## orig 67

	width1 = 67
	##width1=157

rrvar1<- focal(rrvar00, w = matrix(1,width1,width1), fun = fill.na,

           pad = TRUE, na.rm = FALSE)
           
    }       

if(metoden1==0) { rrvar0<-rrvar00 rrvar1<-rrvar0 }

1. 1. summary(getValues(rrvar1))

plot(rrvar1, col=rainbow(8))

  print("invarfname:")
  print(invarfname360_1)
  
  print("Write hadcm3 file ..")

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

writeRaster(rrvar1, invarfname360_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

   rrvar1_180<-rotate(rrvar1)

   rrvar1_180@extent<-extent(-180, 180, -90, 90)

plot(rrvar1_180, col=rainbow(120))

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

writeRaster(rrvar1_180, invarfname180_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

  print("HadCM3 data file loaded .")	

}

1. TPI for different neighborhood size:

tpiw <- function(x, w=15) { m <- matrix(1/(w^2-1), nc=w, nr=w) m[ceiling(0.5 * length(m))] <- 0 f <- focal(x, m) x - f }

1. 1. generate rain shadows, assume westerly wind by default

generate_rain_shadows<-function() { ## params

hillshade_slope=2 hillshade_aspect=260

## direction of global, assumed wind tirektion=(260*6.28)/360

orography1<-raster(outorofname1)

## copy raster windir<-orography1 windir<-windir*0

windir<-windir+tirektion

latras <- lonras <- orography1 xy <- coordinates(orography1) lonras[] <- xy[, 1] latras[] <- xy[, 2]

#print("write windir ...")

crs(windir) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(windir, filename=outname13, overwrite=TRUE, format="CDF", varname="Windir", varunit="radins", longname="Windir", xname="lon", yname="lat")

#print(" write lonras") crs(lonras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(lonras, filename=outname4, overwrite=TRUE, format="CDF", varname="Lonx", varunit="mm", longname="Lonx", xname="lon", yname="lat")

crs(latras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(latras, filename=outname5, overwrite=TRUE, format="CDF", varname="Laty", varunit="mm", longname="Laty", xname="lon", yname="lat")

slope <- terrain(orography1, opt='slope') aspect <- terrain(orography1, opt='aspect') hill <- hillShade(slope, aspect, hillshade_slope, hillshade_aspect) #plot(hill, col=grey(0:100/100), legend=FALSE, main='France')

tpi5 <- tpiw(orography1, w=15)

crs(slope) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(slope, filename=outname7, overwrite=TRUE, format="CDF", varname="Slope", varunit="mm", longname="Slope", xname="lon", yname="lat")

crs(aspect) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(aspect, filename=outname8, overwrite=TRUE, format="CDF", varname="Aspect", varunit="mm", longname="Aspect", xname="lon", yname="lat")

crs(hill) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(hill, filename=outname9, overwrite=TRUE, format="CDF", varname="Hillshade", varunit="m", longname="Hillshade", xname="lon", yname="lat")

crs(tpi5) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(tpi5, filename=outname10, overwrite=TRUE, format="CDF", varname="tp15", varunit="m", longname="tpi5", xname="lon", yname="lat")

#etopo=rout3 #etopo[etopo < 1] <- 0 #etopo[is.na(etopo[])] <- 0

blurelev1 <- focal(orography1, w=matrix(1,nrow=7,ncol=7), fun=mean, pad=TRUE, na.rm = TRUE)

crs(blurelev1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(blurelev1, filename=outname12, overwrite=TRUE, format="CDF", varname="Blurelev", varunit="m", longname="Blurelev", xname="lon", yname="lat")

#plot(rasa1) if(calculate_distance==1) { cacheras1<-orography1

cacheras1[cacheras1 > 1] <- NA cacheras1[cacheras1 < 0] <- 1

print( "Calculating distance. Maybe very slooow. \n Wait ...") distrasa1 <- distance(cacheras1) crs(distrasa1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(distrasa1, filename=outname6, overwrite=TRUE, format="CDF", varname="Dist", varunit="m", longname="Distance", xname="lon", yname="lat") }

}

downscale_dissever <- function (coarse_rastera, fine_stack, dismethod, samplerate) {

   print ("Dissever()")		
   names(fine_stack)
       
     	

coarse_raster<-coarse_rastera

   p1<-fine_stack$Elevation

1. plot(p1)

1. return(0)

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

baset1 <- resample(coarse_raster, p1)

raster_stack <- fine_stack

min_iter <- 3 # Minimum number of iterations max_iter <- 7 # Maximum number of iterations p_train <- samplerate # Subsampling of the initial data

    print("Dissever() begin ...")

oma_juttu <- dissever(coarse = coarse_raster, fine = raster_stack, method = dismethod, p = p_train, min_iter = min_iter,max_iter = max_iter, verbose=1) print("Dissever() end.")

orotemp<-oma_juttu$map

#tempiso<-baset1+oma_juttu$map+biassi

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

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

1. plot(outtemp, col=rev(rainbow(256)) )

1. outtempr<-rotate(outtemp)

#plot(outtempr)

     return(outtemp)

}

downscale_rainfall<-function() { print("Downscaling rainfall ...")

#outname3<-"./data_processing/etopo1.nc" #outname4<-"./data_processing/lons.nc" #outname5<-"./data_processing/lats.nc" #outname6<-"./data_processing/distance.nc" #outname7<-"./data_processing/slope.nc" #outname8<-"./data_processing/aspect.nc" #outname9<-"./data_processing/hillshade.nc" #outname10<-"./data_processing/tpi.nc" #outnam44e11<-"./data_processing/westness.nc" #outname12<-"./data_processing/blurelev.nc" #outname13<-"./data_processing/windir.nc" #outname14<-"./data_processing/etopo_big.nc"

etopo0<<-raster(outorofname1)

  	aspect0<<-raster(outname8)

tpi0<<-raster(outname10) lonx0<<-raster(outname4) hillshade0<<-raster(outname9) slope0<<-raster(outname7) windir0<<-raster(outname13) bluretopo0<<-raster(outname12) latx0<<-raster(outname5)

#print("Krop smallrain")

ext1<-extent(etopo0)

#print(ext1)

#smallrain00<-raster(invarfname180_1) smallrain00<-raster(invarfname180_2)

smallrain0<-crop(smallrain00, ext1)

   plot(smallrain0, col=rev(viridis(100)))	

   aspect1<<-cos(aspect0)

names(aspect1)<<-"Aspect_Cos"

etopo1<<-etopo0 etopo1[etopo1 < 1] <- 1 etopo1<<-log(etopo1) names(etopo1)<<-"Elevation"

landmask0<<-etopo0 landmask0[landmask0 > 0] <<- 1 landmask0[landmask0 < 1] <<- NA

   plot(landmask0)

if(calculate_distance==1) { distans1<<-distans0 distans1[distans1 < 1] <<- 1 distans1<<-log(1/distans1) names(distans1)<<-"Distance" }

lonx1<<-lonx0 lonx1[lonx1 < 1] <<- 1 lonx1<<-log(1/lonx1) names(lonx1)<<-"Lonx"

bhillshade_slope=6 bhillshade_aspect=260

#aspect1=cos(aspect0-windir0)

bslope <<- terrain(bluretopo0, opt='slope') baspect <<- terrain(bluretopo0, opt='aspect')

bhill <<- hillShade(bslope, baspect, bhillshade_slope, bhillshade_aspect)

## windrose w/zulu basis bhill2 <<- ( hillShade(bslope, baspect, bhillshade_slope, bhillshade_aspect)+ hillShade(bslope, baspect, bhillshade_slope, (bhillshade_aspect-20))+ hillShade(bslope, baspect, bhillshade_slope, (bhillshade_aspect+20))+ hillShade(bslope, baspect, bhillshade_slope*2, (bhillshade_aspect-40))+ hillShade(bslope, baspect, bhillshade_slope*2, (bhillshade_aspect+40))+ hillShade(bslope, baspect, 85, (bhillshade_aspect+180)) )/6

bhillshade1<<-bhill halli1<<-(bhillshade1-minValue(bhillshade1))/(maxValue(bhillshade1)-minValue(bhillshade1)) ## zulu basis rainfall correction

#halli2<-(halli1*0.2)+0.9 #halli2=halli1*0.01

1. halli2<<-2+log((0.1+halli1)*10)

   halli2<<-halli1*0.5+0.75

#plot(halli2) names(halli2)<<-"HillshadeX" halli22<<-halli2 names(halli2)<<-"HillshadeX2"

#halli3=log(1+halli2*100) halli3<<-halli2*-1-2

names(halli3)<<-"HSY"

#baspect2=cos(baspect-4.55)*0.05+1.0 #plot(baspect2)

baspect3a<<-((cos(baspect-4.6)+cos(baspect-5.3)+cos(baspect-3.5))/3) baspect2<<-(baspect3a*0.5)+1.0

 #  plot(halli2, col=rev(viridis(100)))
   
 #  plot(baspect2, col=rev(parula(20)))

bhak1<<-baspect2*halli2

   ehak1<<-etopo0*halli2
   ebk1<<-etopo0*baspect2
   ehk1<<-etopo0*bhill
   esk1<<-etopo0*bslope      
   ebhk1<<-etopo0*bhill*baspect2  ##
   ebshk1<<-etopo0*bhill*baspect2*bslope  ## no hyvä

 #  plot(bhak1, col=rev(viridis(100)))
 #  plot(ehak1, col=rev(viridis(100)))
 #  plot(ebk1, col=rev(viridis(100)))
 #  plot(ehk1, col=rev(viridis(100)))
 #  plot(esk1, col=rev(viridis(100)))
 #  plot(ebhk1, col=rev(parula(100))) ## hyvä?
 #        plot(ebshk1, col=rev(viridis(100))) ## no hyvä

names(baspect2)<<-"AspectX" names(ebhk1)<<-"HillShadeAspectX"

#hill <- hillShade(blurslope, bluraspect, hillshade_slope, hillshade_aspect)

tpix1 <<- tpiw(etopo0, w=5)

names(tpix1)<<-"TPIX5"

tpix2 <<- tpiw(etopo0, w=11)

names(tpix2)<<-"TPIX11"

tpix3 <<- tpiw(etopo0, w=15)

names(tpix3)<<-"TPIX15"

tpix4 <<- tpiw(etopo0, w=21)

names(tpix4)<<-"TPIX21"

1. plot(tpix4)

tpix5 <<- tpiw(etopo0, w=51)

names(tpix5)<<-"TPIX51"

plot(tpix5, col=inferno(50) )

tri1<<-terrain(etopo0, opt="TRI")

names(tri1)<<-"TPI1"

plot(tri1, col=viridis(100))

    roughness1<<-terrain(etopo0, opt="roughness")

names(roughness1)<<-"ROUGHNESS1"

plot(roughness1)

    flowdir1<<-terrain(etopo0, opt="flowdir")

names(flowdir1)<<-"FLOWDIR1"

plot(flowdir1)

  dsobject=1

if(dsobject==1) { print("DS 1")

     # preprocess_rasters=1

       plot(etopo0)

g0<-etopo1 g1<-etopo1

#g1<-etopo0*tpix1 #g2<-etopo0*tpix2 #g3<-etopo0*tpix3 #g4<-etopo0*baspect2 g5<-etopo0*halli2 g6<-sin(aspect0)

names(g0)<-"Elevation" names(g1)<-"ElevationX1" #names(g2)<-"ElevationX2" #names(g3)<-"ElevationX3" #names(g4)<-"ElevationX4" #names(g5)<-"ElevationX5" #names(g6)<-"ElevationX6"

       ##rastafari1<-stack(etopo0, icemask_big0, ebhk1) ## kohtal
       # rastafari1<-stack(etopo0, ebhk1) ## kohtal
        rastafari1<-stack(g0,g1) ## kohtal
        

print(names(rastafari1))

#plot(rastafari1$Elevation)

#out3a<<-downscale_dissever(smallrain0, rastafari1,"lm",1.0)

out3<<-downscaler(smallrain0,g0,0)

out3[out3 < 1] <- 1

plot(out3, col=rev(viridis(100)))

print(minValue(out3)) print(maxValue(out3))

#writeout(out3,outvarfname1,"Rainfall", "mm/yr", "IceAgeRain")

       rainout1<-"./data_result/outrain1.nc"

#writeout(out3,rainout1,"Rainfall", "mm/yr", "IceAgeRain")

# out_text_1="Rainfall" # out_text_2="mm/yr" # out_text_3="IceAgeRain"

writeout(out3,rainout1,out_text_1, out_text_2, out_text_3)

  }

}

get_dordogne_srtm<-function() { ## France/Dorgogne area SRTM netcdf file creation ## "R" script

srtm1 <- getData('SRTM', lon=-5, lat=40) srtm2 <- getData('SRTM', lon=-5, lat=45) srtm3 <- getData('SRTM', lon=-5, lat=50) srtm4 <- getData('SRTM', lon=-5, lat=55)

srtm6 <- getData('SRTM', lon=0, lat=40) srtm7 <- getData('SRTM', lon=0, lat=45) srtm8 <- getData('SRTM', lon=0, lat=50) srtm9 <- getData('SRTM', lon=0, lat=55)

srtm10 <- getData('SRTM', lon=5, lat=40)

srtm11 <- getData('SRTM', lon=5, lat=45) srtm12 <- getData('SRTM', lon=5, lat=50) srtm13 <- getData('SRTM', lon=5, lat=55)

srtm14 <- getData('SRTM', lon=10, lat=40) srtm15 <- getData('SRTM', lon=10, lat=45) srtm16 <- getData('SRTM', lon=10, lat=50) srtm17<- getData('SRTM', lon=10, lat=55)

mosaik1 <- merge( srtm1, srtm2,srtm3, srtm4, srtm6,srtm7,srtm8, srtm9,srtm10, srtm11,srtm12,srtm13, srtm14, srtm15,srtm16,srtm17)

#dordogne ext1 <- extent(-2,4,42 , 48)

kropped1<-crop(mosaik1, ext1)

#plot(kropped1)

crs(kropped1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(kropped1, filename='./predata/dordogne.nc', format="CDF", overwrite=TRUE)

   system("del srtm*.*")

}

hadcm3_loadslice <- function(temp_name, var_name, hadyear) {

	putin1 <- nc_open(temp_name)
    	  		          

lones1<<- ncvar_get(putin1, "longitude") latis1<<- ncvar_get(putin1, "latitude") t <- ncvar_get(putin1, "time") lenlones1<-length(lones1) lenlatis1<-length(latis1)

deltayears1=hadyear-hadbaseyear deltamonths1=deltayears1*12 item1=30000-deltamonths1-12+1 months1=12

temp_pusu1<-ncvar_get(putin1,var_name, start=c(1,1,item1), count=c(lenlones1,lenlatis1,months1) ) nc_close(putin1)

taimi1=t[1]

return(temp_pusu1)

}

generate_hadfilename<-function(hadbaspath1, yrr1, varr1) {

hadfilenamex1=hadbaspath1

hadfilenamex1<-paste0(hadfilenamex1,"/bias_regrid_") hadfilenamex1<-paste0(hadfilenamex1,varr1) hadfilenamex1<-paste0(hadfilenamex1,"_")

a=as.integer(yrr1/2500) b=a*2500 c=b/1000 d=c+2.5

hadbaseyear<<-b

hadfilenamex1<-paste0(hadfilenamex1,toString(c)) hadfilenamex1<-paste0(hadfilenamex1,"_") hadfilenamex1<-paste0(hadfilenamex1,toString(d)) hadfilenamex1<-paste0(hadfilenamex1,"kyr.nc") return(hadfilenamex1) }

load_had_slices<-function(beginyr1, yrs1, varr1) { endyr1=beginyr1+yrs1-1 print("Loading haccm3 slices, wait ...")

markki1=0 yyyy1=0

for (yrr1 in (beginyr1:endyr1)) { hadfilename=generate_hadfilename(hadbasepath, yrr1, varr1) print(yrr1) print (hadfilename) slice00=hadcm3_loadslice(hadfilename, varr1, yrr1)

if(markki1==0) { baseslice1<-slice00 } else { # add slices baseslice1<-baseslice1+slice00 }

markki1=1 yyyy1=yyyy1+1 }

		#print(head(baseslice1))
		baseslice1=baseslice1/yyyy1
		

return (baseslice1) }

draw_climate_diagram<-function(lampot, sadem) {

#mydata <- read.csv("kiova2.txt", header=FALSE, sep=";") labeli='Paris, 40750 BP' nimi="paris_40750bp" datanimi=paste(nimi,".txt"); kuvanimi=paste(nimi,".svg");

prmax=100 prmin=0 tmax = 20.0 tmin=-25.0 tstep=5

widthi=10 heighti=16

asteikko<-c(" "," ","3"," "," ","6"," ", " ","9"," "," ","12" )

svg(kuvanimi, width=widthi, height=heighti)

deltapr=prmax-prmin deltatee<-(tmax-tmin)

ratio<-deltapr/deltatee

y2offset= -1*ratio*tmin

total_sadem=sum(sadem) avg_lampotila=sum(lampot)/12 avg_lampotila=(round(avg_lampotila)*10)/10

max_lampotila=max(lampot) min_lampotila=min(lampot)

par(mar=c(6,6,6,6),cex.axis=2,cex.lab=2.5)

b<-barplot(sadem, names.arg=asteikko, col="blue", border="blue",ylim=c(prmin, prmax), cex.axis=2.5, cex.names=2.5 )

lines(b, (lampot*ratio)+y2offset, col="Red",lwd=8)

right.axis.ticks<- seq(from =tmin , to=tmax , by=tstep)

axis(4,at=(right.axis.ticks*ratio)+y2offset,labels=paste0(right.axis.ticks),las=2, cex.axis=2.5)

mtext(side = 2, line = 3, 'Precipitation', cex=2.5, col="darkblue") mtext(side = 4, line = 3, 'Temperature', cex=2.5, col="darkred") mtext(side = 1, line = 3, 'Month', cex=2.5, col="darkgreen")

text(7,(prmax-2),cex=3.5, labeli); text(1,(prmax-8),cex=2.4, pos=4, paste("Tavg=",avg_lampotila, " C" )); text(1,(prmax-12),cex=2.4, pos=4, paste("Tamax=",max_lampotila, " C" )); text(1,(prmax-16),cex=2.4, pos=4,paste("Tamin=",min_lampotila, " C" )); text(1,(prmax-20),cex=2.4, pos=4,paste("Pra=",total_sadem, " mm" ));

}

get_had_climate_data<-function(beginyears, years, targetname1, lat1, lon1) { print("Loading data, wait ...") varr1="tas" varr2="pr" tempsit1<-load_had_trapezoid(beginyears, years, varr1, lon1, lat1) precsit1<-load_had_trapezoid(beginyears, years, varr2, lon1, lat1)

#print (tempsit1) #print (precsit1)

months1<-1:12

tempsit1<-round(tempsit1, digits = 1) precsit1<-round(precsit1, digits = 1)

tavg1<-sum(tempsit1)/12.0 pannual1<-sum(precsit1)

df1<-data.frame(months1, tempsit1,precsit1)

names(df1)<-c("Month", "T", "P")

coutname1=paste0(targetname1, ".csv")

#write.csv2(df1,coutname1)

#write.table(df1,file=coutname1,sep=";")

   write.table(df1,file=coutname1,sep=";",row.names=FALSE)

print("Monthly data:") print(df1)

print ("Climate averages:") print (tavg1) print (pannual1)

draw_climate_diagram(tempsit1, precsit1)

}

had_twoslicer<-function(beyr1,yrs1,varr1) {

enyr1=beyr1+yrs1

print(beyr1) print(enyr1)

hadbasepath1<<-hadbasepath

hadnames1<-vector(mode="character", length=2)

hadbaseyears<-rep(0,2) #print (hadbaseyears)

hadnames1[1]<-generate_hadfilename(hadbasepath1, beyr1, varr1) hadbaseyears[1]<-hadbaseyear hadnames1[2]<-generate_hadfilename(hadbasepath1, enyr1, varr1) hadbaseyears[2]<-hadbaseyear

deltayears1=(beyr1-hadbaseyears[2])*-1 deltamonths1=deltayears1*12 item1=30000-deltamonths1-12+1 months1=12

deltayears2=enyr1-hadbaseyears[2] deltamonths2=deltayears2*12 item2=30000-deltamonths2-12+1 months2=12

1. print (hadnames1[1])
2. print (hadnames1[2])
3. print(deltayears1)
4. print(deltamonths1)
5. print(deltayears2)
6. print(deltamonths2)

twoo1=0

if(hadbaseyears[1]==hadbaseyears[2]) {

1. print("Twoo 1")

twoo1=1 }

if(deltayears1>-1) { twoo1=1 }

if(twoo1==0) {

putin1 <- nc_open(hadnames1[1])

lones1<<- ncvar_get(putin1, "longitude") latis1<<- ncvar_get(putin1, "latitude") t <- ncvar_get(putin1, "time") lenlones1<-length(lones1) lenlatis1<-length(latis1)

temp_pusu1<-ncvar_get(putin1,varr1, start=c(1,1,item1), count=c(lenlones1,lenlatis1,deltamonths1) ) nc_close(putin1)

1. print("put in 2")

putin2 <- nc_open(hadnames1[2])

lones1<<- ncvar_get(putin2, "longitude") latis1<<- ncvar_get(putin2, "latitude") t2 <- ncvar_get(putin2, "time") lenlones1<-length(lones1) lenlatis1<-length(latis1)

temp_pusu2<-ncvar_get(putin2,varr1, start=c(1,1,item2), count=c(lenlones1,lenlatis1,deltamonths2) ) nc_close(putin2)

# print("put in 2")

dima1=dim(temp_pusu2)

pusu3=abind(temp_pusu1,temp_pusu2,along=3)

} #two file buffers else { # print("Twoo 1 ...")

deltayears2=beyr1-hadbaseyears[1] deltamonths2=deltayears2*12 item2=30000-deltamonths2-12+1 months2=(enyr1-beyr1)*12

#print(deltayears2) #print(deltamonths2) #print(item2)

#print("put in 2") putin2 <- nc_open(hadnames1[2])

lones1<<- ncvar_get(putin2, "longitude") latis1<<- ncvar_get(putin2, "latitude") t2 <- ncvar_get(putin2, "time") lenlones1<-length(lones1) lenlatis1<-length(latis1)

pusu3<-ncvar_get(putin2,varr1, start=c(1,1,item2), count=c(lenlones1,lenlatis1,months2) ) nc_close(putin2)

}

dima3=dim(pusu3)

#print(dima1) #print(dima3)

as1<- array(rep(0, 720*180*12), dim=c(720, 180, 12))

ylimit1=dima3[3]

#print (dim(as1))

hhh1=0

maxima1<-(ylimit1/12)-1

print (maxima1) for( m in 1:maxima1) { for( n in 1:12) {

has1<-pusu3[,,m*12+n]

as1[,,n]<-as1[,,n]+pusu3[,,m*12+n]

} hhh1=hhh1+1 }

as1<-as1/hhh1

return(as1) }

load_had_trapezoid<-function(beginyr1, yrs1, varr1, lon1, lat1) {

1. slice00=load_had_slices(beginyr1, yrs1, varr1)

slice00<-had_twoslicer(beginyr1,yrs1,varr1)

dima1=dim(slice00)

#print (dima1)

max1=dima1[1] may1=dima1[2]

   londex2=which(lones1 >= lon1 )[1]
   latdex2=which(latis1 >= lat1 )[1]
  
  
   londex1=londex2-1
   latdex1=latdex2-1
  
   if(londex1<1) londex1=max1	
   if(latdex1<1) latdex1=may1
     	
   abslon1=lones1[londex1]
   abslat1=latis1[latdex1]	
   abslon2=lones1[londex2]
   abslat2=latis1[latdex2]	

#print("lons") #print(abslon1)

   #print(abslon2)
   #print(abslat1)
   #print(abslat2)
   
   #print (max1)
   #print (may1)
   #print (lones1[0])	

vektor1<-1:12 vektor1<-vektor1*0

   n=7
  	for (n in 1:12)

{ ## attempt to process trapezoid

value1=slice00[londex1,latdex1, n] value2=slice00[londex1,latdex2, n] value3=slice00[londex2,latdex1, n] value4=slice00[londex2,latdex2, n]

rulon1=abslon2-abslon1 rulat1=abslat2-abslat1

daata1<-c(value1,value2,value3,value4)

matrix <- matrix(daata1, nrow=2, ncol=2) r <- raster(matrix) ## lon lat extent(r) <- c(abslon1, abslon2, abslat1,abslat2)

## reso 100x100 s <- raster(nrow=100, ncol=100)

extent(s)<-extent(r) s <- resample(r, s, method='bilinear')

xy <- cbind(lon1,lat1)

resultt1<-extract(r, xy)

vektor1[n]=resultt1 }

return(vektor1)

}

load_had_raster<-function(beginyr1, yrs1, varr1, month1) { #slaici1=load_had_slices(beginyear1, yrs1, varr1) slaici1<-had_twoslicer(beginyr1,yrs1,varr1)

dima1=dim(slaici1)

print (dima1)

if(month1==0) { markki=0 yyyy1=0

## select all months for (n in 1:12) { vaari0=slaici1[,,month1] if(markki1==0) { baseslice1<-slice00 } else { # add slices baseslice1<-baseslice1+slice00 }

merkki=1 yyyy1=yyyy1+1 }

vaari0=baseslice1/yyyy1 } else { vaari0=slaici1[,,month1] }

print (dim(vaari0))

padding1 = matrix(0,720,180)

   vaari1<-cbind(padding1,vaari0)  

vaari1<- apply(t(vaari1),2,rev)

rrvar1<-raster (vaari1)

rrvar1@extent<-extent(0, 360, -90, 90)

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

   rvarfilename1=paste0(hadprocesspath, "global_360_",varr1,"_",month1,"-nc")
   longname1=paste0(varr1," ",toString(beginyr1) )

writeRaster(rrvar1, rvarfilename1, overwrite=TRUE, format="CDF", varname=varr1, varunit="unit", longname=longname1, xname="lon", yname="lat")

}

load_had_rasters_var<-function(byr1, yrs1, varr1) { yrmid1=byr1+(yrs1/2)

slaici1<-had_twoslicer(byr1,yrs1,varr1)

dima1=dim(slaici1)

print (dima1)

for (n in 1:12) { print (n) vaari0=slaici1[,,n]

padding1 = matrix(0,720,180)

vaari1<-cbind(padding1,vaari0)

vaari1<- apply(t(vaari1),2,rev)

rrvar1<-raster (vaari1)

rrvar1@extent<-extent(0, 360, -90, 90)

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

if(n==1) { rs1=stack(rrvar1) } else { rs1=stack(rs1, rrvar1) }

}

 plot(rs1)

   rvarfilename1=paste0(hadprocesspath, "global_360_",varr1,"_",yrmid1)
   longname1=paste0(varr1," ",toString(yrmid1) )

writeRaster(rs1, rvarfilename1, overwrite=TRUE, format="CDF", varname=varr1, varunit="unit", longname=longname1, xname="lon", yname="lat")

   inputdataname1<<-rvarfilename1
    

}

load_climate<-function() {

1. beginyear1=40700

beginyear1=33950

years1=100

month1=7

varr1="tas" varr2="pr"

1. 1. paris
2. targetname1="paris"
3. targetlat1=48.856667
4. targetlon1=2.351111

targetname1="sungir" targetlat1=56.175833 targetlon1=40.509167

1. targetlon1=0.0

print("-----------------------------") print("Age:") hage1<-beginyear1+(years1/2) print(hage1) print("Target:") print(targetname1) print (targetlon1) print (targetlat1) print ("")

get_had_climate_data(beginyear1, years1,targetname1,targetlat1, targetlon1)

}

hadcm3_raster_fillna_method_2<-function(vaari3) { printf(" FillNA method 2. Sooow, wait patiently ...") vaari4<-vaari3 vaari5<-vaari3 vaari6<-vaari3

   	dima1<-dim(vaari3)

w1=dima1[1] h1=dima1[2] #print (w1) #print (h1) for (ix in (1:w1) )

{ beg1=1 viksu1=-1 for (iy in (1:h1) )

{ piksu1=vaari3[ix,iy]

if(is.na(piksu1)) {

if(beg1==0) { vaari4[ix,iy]=viksu1 vaari3[ix,iy]=1 } else { vaari4[ix,iy]=piksu1

}

} else { beg1=0 viksu1=piksu1 beg1=0 }

}

beg1=1 viksu1=-1 for (iy in seq(h1, 1, by=-1 ) ) {

piksu1=vaari3[ix,iy]

if(is.na(piksu1)) {

vaari4[ix,iy]=viksu1

} else ## ei na { beg1=0 viksu1=piksu1 beg1=0 }

} }

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. sekond

  				for (iy in (1:h1) )
   

{ beg1=1 viksu1=-1 for (ix in (1:w1) )

{ piksu1=vaari5[ix,iy]

if(is.na(piksu1)) {

if(beg1==0) { vaari6[ix,iy]=viksu1 vaari5[ix,iy]=1 } else { vaari6[ix,iy]=piksu1

}

} else { beg1=0 viksu1=piksu1 beg1=0 }

}

beg1=1 viksu1=-1 for (ix in seq(w1, 1, by=-1 ) ) {

piksu1=vaari5[ix,iy]

if(is.na(piksu1)) {

vaari6[ix,iy]=viksu1

} else ## ei na { beg1=0 viksu1=piksu1 beg1=0 }

} }

vaari7<-(vaari4+vaari6)/2 return(vaari7) }

fill.na <- function(x) { width1<-67

 center = 0.5 + (width1*width1/2) 
 if( is.na(x)[center] ) {
   return( round(mean(x, na.rm=TRUE),0) )
 } else {
   return( round(x[center],0) )
 }

}

fill.na7 <- function(x) { width1<-7

 center = 0.5 + (width1*width1/2) 
 if( is.na(x)[center] ) {
   return( round(mean(x, na.rm=TRUE),0) )
 } else {
   return( round(x[center],0) )
 }

}

get_hadgcm3_global_raster<-function(byr1, yrs1, var1, month1, varunit1, ovarnam1, ofisuffix1, method1) {

1. 1. mont 0 sum of all year
   2. month 13 average of all year
   3. method 4: resample ratser in input, attempting
   4. toreduce original downscaling out

   print ("Get global raster ...")
   

load_had_rasters_var(byr1, yrs1, var1)

inputdataname10<-paste0(inputdataname1,".nc")

stak1<-stack(inputdataname10)

if(month1==0) { r00 <- calc(stak1, sum) }

if(month1==13) { r00 <- calc(stak1, sum)/12 } else { if(month1>0) { r00<-stak1@layersmonth1 } }

oname1=paste0(hadprocesspath,"/",ofisuffix1,"0_360.nc") writeout(r00, oname1, var1, "C", var1)

rr00_180<-rotate(r00)

## rr00_180@extent<-extent(-180, 180, -90, 90)

plot(rr00_180, col=rainbow(64))

oname2=paste0(hadprocesspath,"/",ofisuffix1,"0_180.nc") writeout(rr00_180, oname2, var1, "C", var1)

inras1<-raster(oname1) inras2<-raster(oname2) #varunit1="C" if(method1==0) {

sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc") writeout(inras1,sampledname1, var1, varunit1, var1) sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc") writeout(inras2,sampledname2, var1, varunit1, var1) }

if(method1==1) { ## orig 67 width1 = 67 ##width1=157

outras1<- focal(inras1, w = matrix(1,width1,width1), fun = fill.na, pad = TRUE, na.rm = FALSE)

       outras2<- focal(inras2, w = matrix(1,width1,width1), fun = fill.na, pad = TRUE, na.rm = FALSE)   

sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc") writeout(outras1,sampledname1, var1,varunit1, var1) sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc") writeout(outras2,sampledname2, var1,varunit1, var1) }

  if(method1==2)
  { 

outras1<-hadcm3_raster_fillna_method_2(inras1) outras2<-hadcm3_raster_fillna_method_2(inras2) sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc") writeout(outras1,sampledname1, var1,varunit1, var1) sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc") writeout(outras2,sampledname2, var1,varunit1, var1) }

if(method1==4) { sampler1 <- raster(ncol=samplecols1, nrow=samplerows1) outras01<-resample(inras1, sampler1) outras02<-resample(inras2, sampler1)

width1=7

outras1<- focal(outras01, w = matrix(1,width1,width1), fun = fill.na7, pad = TRUE, na.rm = FALSE)

       outras2<- focal(outras02, w = matrix(1,width1,width1), fun = fill.na7, pad = TRUE, na.rm = FALSE) 

plot(outras1)

sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc") writeout(outras1,sampledname1, var1,varunit1, var1) sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc") writeout(outras2,sampledname2, var1,varunit1, var1) }

}

test_code_raster_experiment_1<-function() { beginyear1=45000-5 years1=10 month1=7 varr1="tas" varr2="pr" method1=4 ovarnam1="TS" varunit1="C" ofisuffix1="tempin"

#get_hadgcm3_global_raster(beginyear1, years1, varr1, month1, method1)

get_hadgcm3_global_raster(beginyear1, years1, varr1, month1,varunit1, ovarnam1,ofisuffix1,method1) }

1. 2. Main proggis

1. 1. main running ...

print("TS downscaler.")

print ("Main program running ...")

read_argus()

print("Age") print(fage) print("Numyears") print(numyears) print("Months") print(fmonth)

sealevel_from_age(fage)

if (orography_preprocess==1) { "Orography ..." preprocess_orography() }

if(get_srtm_data==1) { get_dordogne_srtm() }

if(get_gtopo30_data==1) { get_europe_gtopo30() }

if(capture_hadcm3_60ka_data==1) {

if(capture_temperature==1)

{
capture_hadcm3_temperature=1
}

if(capture_rainfall==1)
{
capture_hadcm3_rainfall=1
}

if(capture_hadcm3_temperature==1)

{ print("Processing hadcm3 temperature data to snapshot ...")

hadvariaabeli="tas"

	#hadfilename<<-generate_hadfilename(hadbasepath, hadyear, hadvariaabeli)

#print (hadyear)

#load_average_hadcm3_years_months(hadfilename, hadvariaabeli,hadbaseyear, hadyear,hadmonth, hadyears, hadmonths, hadoperaatio,hadcm3_60ka_seafill_method)

#byr1=hadbaseyear-(hadyears/2) byr1=fage years1=numyears month1=fmonth varr1=hadvariaabeli method1=hadcm3_60ka_seafill_method ovarnam1="TS" varunit1="C" ofisuffix1="tempin"

#print ("BYR") #print (byr1) get_hadgcm3_global_raster(byr1, years1, varr1, month1,varunit1, ovarnam1,ofisuffix1,method1)

}

if(capture_hadcm3_rainfall==1)

{

print("Processing hadcm3 rainfall data to snapshot ...") hadvariaabeli="pr"

#hadfilename<<-generate_hadfilename(hadbasepath, hadyear, hadvariaabeli)

#load_average_hadcm3_years_months(hadfilename, hadvariaabeli,hadbaseyear, hadyear,1, hadyears, 12, 3, hadcm3_60ka_seafill_method)

byr1=fage years1=numyears month1=fmonth varr1=hadvariaabeli method1=hadcm3_60ka_seafill_method ovarnam1="PR" varunit1="C" ofisuffix1="precipin"

#print ("BYR") #print (byr1) get_hadgcm3_global_raster(byr1, years1, varr1, month1,varunit1, ovarnam1,ofisuffix1,method1)

  }

}

if(capture_trace26k_data==1) { if(capture_temperature==1) { print("Processing input temperature variable to snapshot ...") tracelocation3(predata2, variaabeli, fage, numyears, fmonth); } if(capture_rainfall==1) { print("Processing input rainfall variable to snapshot ...") tracelocation3(predata2, variaabeli, fage, numyears, fmonth); } }

if(capture_elevation==1) { print("Glac elevation ...") load_glac_elev() load_glac_icemask() }

if(capture_elevation==2) { print("Peltier elevation") load_peltier_elev() }

if(make_rain_shadows==1) { if(capture_elevation==2)

  {	

# force method to 0 because of peltier topography downscale_orography(0)

  }
  else
  {	

downscale_orography(method_ds_oro)

  }

## generate rain shadows, assume westerly wind by default print("Generating rain shadows.") generate_rain_shadows()

}

print ("Downscaling ...")

1. 1. do ts downscaling

if(normal_ds==1) { print("Normal etopo1 area downscaling.") ## default area DS normal_ts_downscale() }

if(normal_ds==3) { print("Rainfall downscaling.")

#out_text_1="Rainfall" #out_text_2="mm/yr" #out_text_3="IceAgeRain"

## default area DS downscale_rainfall() }

if(global_ds==1) { print("World downscaling.") world_ts_downscale() }

if(accurate_ds==1) { print ("Customized dwonscaling.") custom_ts_downscale_1(kustomorofilee1,outvarfname1 ) }

if(plot_result==1) { print("Plotting of draft.")

 plottaus1(outvarfname1,plotfname1 )

}

1. try(system("python maplot1.py", intern = TRUE, ignore.stderr = TRUE))

print (fage) print (fmonth) print("Program run done.")