date<-strptime(sampling_dat$date[i], "%Y%m%d")   # interpolation date being processed, converting the string using specific format

  LST_month<-paste("mm_",month,sep="") # name of LST month to be matched in the raster stack of covariates and data.frame

  pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA

  s_raster<-dropLayer(s_raster,pos)      # If it exists drop layer

  ghcn.subsets[[i]] <- transform(ghcn.subsets[[i]],LST = mod_LST)            #Add the variable LST to the subset dataset

  dst$LST<-dst[[LST_month]] #Add also to monthly dataset

  #  STEP 1 - LST 10 year monthly averages: THIS IS NOT USED IN CAI method

  names(modst)[1]<-c("id")       #modst contains the average tmax per month for every stations...it has 193 rows

  dmoday=merge(modst,d,by="id",suffixes=c("",".y2"))  #LOOSING DATA HERE!!! from 113 t0 103

  xmoday=merge(modst,x,by="id",suffixes=c("",".y2"))  #LOOSING DATA HERE!!! from 48 t0 43

  mod_pat<-glob2rx("*.y2")   

  var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names

  dmoday<-dmoday[,-var_pat]

  mod_pat<-glob2rx("*.y2")   

  var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names

  xmoday<-xmoday[,-var_pat] #Removing duplicate columns

  #dmoday=merge(modst,d,by="id")  #LOOSING DATA HERE!!! from 113 t0 103

  #xmoday=merge(modst,x,by="id")  #LOOSING DATA HERE!!! from 48 t0 43

  #names(dmoday)[4]<-c("lat")

  #names(dmoday)[5]<-c("lon")     #dmoday contains all the the information: BIAS, monn

  #names(xmoday)[4]<-c("lat")

  #names(xmoday)[5]<-c("lon")     #dmoday contains all the the information: BIAS, monn

  clim_xy<-project(as.matrix(sta_lola),proj_str)     #This is the coordinates of monthly station location (193)

    sta_bias<-dmoday$LSTD_bias         

    bias_xy<-cbind(dmoday$x_OR83M,dmoday$y_OR83M) #This will use only stations from training daily samples for climatology step if bias_val=1

  sta_clim<-modst$TMax #This contains the monthly climatology...used in the prediction of the monthly surface

  # STEP 7 - interpolate delta across space: this is the daily deviation from the monthly average

  daily_deltaclim<-dmoday$dailyTmax-dmoday$TMax    #For daily surface interpolation...

  daily_deltaclim_v<-data_v$dailyTmax-data_v$TMax  #For validation...

  # STEP 8 - assemble final answer - T= LST-Bias(interpolated)+delta(interpolated)    (This is for fusion not implemented in this script...)

  #                                  T= clim(interpolated) + deltaclim(interpolated)  (This is for CAI)

  coords<- modst[,c('x_OR83M','y_OR83M')]

  coordinates(modst)<-coords

  proj4string(modst)<-CRS  #Need to assign coordinates..

  if (climgam==1){          #This is an option to use covariates in the daily surface...

    data_month<-modst

    data_month$y_var<-modst$TMax

  #mod1<- try(gam(formula1, data=data_s))

  #mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2

  #mod3<- try(gam(formula3, data=data_s))

  #mod4<- try(gam(formula4, data=data_s))

  #mod5<- try(gam(formula5, data=data_s))

  #mod6<- try(gam(formula6, data=data_s))

  #mod7<- try(gam(formula7, data=data_s))

  #mod8<- try(gam(formula8, data=data_s))

  if (climgam==1){          #This will automatically use monthly station data in the second step

    mod1<- try(gam(formula1, data=data_month))

    mod2<- try(gam(formula2, data=data_month)) #modified nesting....from 3 to 2

    mod3<- try(gam(formula3, data=data_month))

    mod4<- try(gam(formula4, data=data_month))

    mod5<- try(gam(formula5, data=data_month))

    mod6<- try(gam(formula6, data=data_month))

    mod7<- try(gam(formula7, data=data_month))

    mod8<- try(gam(formula8, data=data_month)) 

  } else if (climgam==0){ #This will use daily delta in the second step

    mod1<- try(gam(formula1, data=data_s))

    mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2

    mod3<- try(gam(formula3, data=data_s))

    mod4<- try(gam(formula4, data=data_s))

    mod5<- try(gam(formula5, data=data_s))

    mod6<- try(gam(formula6, data=data_s))

    mod7<- try(gam(formula7, data=data_s))

    mod8<- try(gam(formula8, data=data_s))

  }

  #browser()

      # model specific metrics

        tmax_predicted_CAI<-raster_pred + clim_rast #Final surface  as a raster layer...taht is if daily prediction with GAM

        raster_name<-paste("GAMCAI_tmax_predicted_",data_name,out_prefix,".rst", sep="")

      #y_var_fit= mod$fit #move it

      #Use res_mod_s so the R2 is based on daily station training

      R2_mod_f<- cor(data_s$dailyTmax,res_mod_s, use="complete")^2

      RMSE_mod_f<- sqrt(mean(res_mod_s^2,na.rm=TRUE))

      results_RMSE_f[1,1]<- sampling_dat$date[i]  #storing the interpolation dates in the first column

      results_RMSE_f[1,2]<- ns        #number of stations used in the training stage

      results_RMSE_f[1,3]<- "RSME_f"

      #results_RMSE_f[1,j+3]<-sqrt(mean(mod$residuals^2,na.rm=TRUE))

      results_RMSE_f[1,j+3]<-sqrt(mean(res_mod_s^2,na.rm=TRUE))

      results_MAE_f[1,1]<- sampling_dat$date[i]  #storing the interpolation dates in the first column

      results_MAE_f[1,2]<- ns        #number of stations used in the training stage

      results_MAE_f[1,3]<- "MAE_f"

      #results_MAE_f[j+3]<-sum(abs(y_var_fit-data_s$y_var))/ns

      results_MAE_f[1,j+3]<-mean(abs(res_mod_s),na.rm=TRUE)

      results_R2_f[1,1]<- sampling_dat$date[i]      #storing the interpolation dates in the first column

      results_R2_f[1,2]<- ns            #number of stations used in the training stage

      results_R2_f[1,3]<- "R2_f"

      results_R2_f[1,j+3]<- R2_mod_f      #Storing R2 for the model j

      #### Now calculate validation metrics

  #results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj)

  results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,data_month)

  names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","data_month")