

#################################################################################.
#                             :: EBC_Pipeline.R ::
# 
# Pipeline for the Energy Balance Closure analysis in ICOS and NEON
# author:          Giacomo Nicolini
# creation date:   20241018
# last update:     20260323
#################################################################################.

proj.dir <- getwd()

options(scipen = 999)


# Libraries ---------------------------------------
kpacks <- c("data.table",
            "maps",
            'RColorBrewer',
            'zoo',
            'utils',
            'stringr',
            'lubridate',
            'dplyr'
)
new.packs <- kpacks[!(kpacks %in% installed.packages()[,"Package"])]
if(length(new.packs)){install.packages(new.packs)} # install packages not already installed
loaded.packs <- lapply(kpacks, require, character.only=T) # load them
if(!all(unlist(loaded.packs))){cat(yellow $ bold(paste('WARNING: not able to load the following libraries:\n', paste(kpacks[which(unlist(loaded.packs)==F)], collapse=', '))))}
remove(kpacks, new.packs)

## Functions ---------------------------------------
# Local functions ..............................
source('Functions/Solar_zenith_time_functions.R') # A set of functions to compute the solar zenith angle and daytime/nighttime/twilight (used to retrieve "solar_zenith" function only)
                                                  # (used to retrieve "solar_zenith" and "classify_daynight" functions only)
source('Functions/SWIN_slope_correction_v3.R')    # A set of functions to compute the radiation topographic corrections
source('Functions/RadOUT_tower_disturbance.R')    # A set of functions to estimate the disturbance of the tower structure on the outgoing SW and LW radiations and visualizations
source('Functions/Make_combos_filtered.R')        # Function to make and filter AE combos
source('Functions/EBC_calc.R')                    # Function to perform regression among AE and TE

## Processing settings -------------
# Set the footprint cumulative contribution threshold to filter out data in the TE6 scenario
ffp.thrhld <- 0.6

# Set the conversion factors for the conversion of GPP from μmolCO2m-2s-1 to kgCO2m-2s-1 (used for photosynthetic energy consumption)
umol2mol <- 10^-6   
CO2.mm <- 0.044     # molar mass of CO2
pecf <- -10.88*10^6 # photosynthetic energy conversion factor  J kgCO2-1

# Set the canopy-specific heat capacity = 2958 J K−1 kg−1 (used for heat stored in biomass)
cshc <- 2958 # 2000 in Kuticoff

# Set the "net" radiometer FoV angle
rad.FoV_SWOUT <- 170
rad.FoV_LWOUT <- 150

# Set the albedo of the tower (see Kidston et al. 2010)
# Acciaio grezzo / metallo ossidato (ferro, acciaio non verniciato, ruggine): Albedo 0.08 – 0.25
# Acciaio verniciato (chiaro / bianco): Albedo 0.6 – 0.9
albedo_tower <- 0.2

# Set the emissivity of the tower (see Kidston et al. 2010)
emissivity_tower <- 0.3

# New dataset: Spho, Sbio, NR corrected data, ...
save.newdata <- TRUE

# Plotting controls
save.RAD_slope_correction <- TRUE       # whether to save the SW_IN slope correction DC plot
save.RAD_dem_correction <- TRUE         # whether to save the topographic shadowing correction DC plot
save.RAD_towdist_omega <- TRUE          # whether to save the RadFoV-tower correction solid angles scheme 
save.RAD_towdist_correction <- TRUE     # whether to save the RadFoV-tower correction plots 
save.RAD_towdist_emissivity <- TRUE     # whether to save the RadFoV-tower correction estimated emissivity vs LW_OUT
save.RAD_fovffp_correction <- TRUE      # whether to save the Rad vs FFP FoV correction plots 
save.NETRAD_full_correction <- FALSE     # whether to save the daily cycle of NETRAD after complete correction (all of the above)  
save.spcomp <- FALSE                     # whether to save the spectral correction comparison


# processing version number
vernum <- '16'



# Ancillary data import --------------------------

# . STATIC ANCILLARY data table (with ecosystem, biomass, coordinates, slope, profile, etc.)
ebc.ancillary <- data.frame()
ebc.ancillary <- fread(paste0(proj.dir, '/EBC_DATA/EBC_ancillary_data/EBC_stations_MD_ANCILLARY.txt'), sep='\t', header = T, data.table = F)

# . DYNAMIC ANCILLARY data lists (biomass and LAI by station, (LAI ICOS only))
bio.data.dyn.icos <- readRDS("EBC_DATA/EBC_ancillary_data/EBC_biomass_dynamic_ICOS.RData")
bio.data.dyn.neon <- readRDS("EBC_DATA/EBC_ancillary_data/EBC_biomass_dynamic_NEON.RData")

# . DYNAMIC canopy heights data lists (by station, ICOS only)
hc.data.dyn <- readRDS("EBC_DATA/EBC_ancillary_data/EBC_canopyheight_dynamic_ICOS.RData")

# . TA layer-weighted
# TA layer-weighted values 
TA.lw.files <- list.files('C:/WORK/R/PROJECTS/EBC-ICOS/EBC_DATA/EBC_ancillary_data/TA_layer-weighted/', full.names = TRUE)

# . Radiation setup 
# table done with the inputs of station PIs
# h_radiometro: altezza radiometro dal suolo
# w_traliccio: lato traliccio portante radiometro
# w_traliccio_net: lato traliccio portante radiometro al netto degli spazi vuoti, considerati al 40%
# d_braccio: lunghezza boom radiometro, sul piano ortogonale al lato del traliccio (corretto se alpha != 90°)
rad.setup <- fread(paste0(proj.dir, '/EBC_DATA/EBC_ancillary_data/ICOS_stations_radiometer_setup.txt'), sep='\t', header = T, data.table = F)
rad.setup[rad.setup == 0] <- NA # put 0 values to NA

# . EC-RAD FOVs 
# Table with the ratios of 'radiations' measured by the radiometer and those measured with EC
# rasters of visible+IR bands have been cropped to the FFP climatology (70%) and to the radiometer FOV at the canopy top
# the indexes were weighed by the respective pixel contributions (footprint values, constant values for radiometer) and integrated over the respective areas.
# These values (RAD/FFP) have been related (r > 1 means the radiometer 'overestimate' with respect to the EC fov (footprint))
# 2024 data only, data are split by months
# P: N progressivo periodo 30 giorni
# vis_nir: rapporto FFP/RADFoV media bande visibile, NIR e SWIR, da usare per correzione SW_OUT
# surf_temp: rapporto FFP/RADFoV media bande per suface temperatur, da usare per correzione LW_OUT
# diff_area: rapporto tra l'intersezione FFP/RADFoV e l'area maggiore tra FFP e RADFoV, per quantificare la differenza fisica ra le due
rad.ffp.fov <- fread(paste0(proj.dir, '/EBC_DATA/EBC_ancillary_data/ICOS_stations_RADvsFFP_fov.csv'), sep=',', header = T, data.table = F)

# Data import --------------------------

# DATA files base dir
data.dir <- "EBC_DATA/EBC_input_data/"

# ICOS
ebc.icos.data <- list()
ebc.icos.data <- readRDS(paste0(data.dir , "EBC_data_ICOS.RData"))
flxnt.icos.data <- list()
flxnt.icos.data <- readRDS(paste0(data.dir , "Fluxnet_data_ICOS.RData"))
flxnt.icos.usttrhd <- readRDS(paste0(data.dir , "Fluxnet_ustar_ICOS.RData"))
ebc.GF.icos.data <- list()
ebc.GF.icos.data <- readRDS(paste0(data.dir , "EBC_data_ICOS_GF.RData"))
nrt.icos.data <- list()
nrt.icos.data <- readRDS(paste0(data.dir , "NRT_data_ICOS.RData"))

# NEON
ebc.neon.data <- list()
ebc.neon.data <- readRDS(paste0(data.dir , "EBC_data_NEON.RData"))
flxnt.neon.data <- list()
flxnt.neon.data <- readRDS(paste0(data.dir , "Fluxnet_data_NEON.RData"))
flxnt.neon.usttrhd <- readRDS(paste0(data.dir , "Fluxnet_ustar_NEON.RData"))
ebc.GF.neon.data <- list()
ebc.GF.neon.data <- readRDS(paste0(data.dir , "EBC_data_NEON_GF.RData"))



#  *** EBC BASE DATA ***  -----------------------------------------
# a list of dataframes containing the timeseries of Available Energy (AE) and Turbulent Energy (TE) scenarios per site

# Trim the station/ancillary table to the effective data availability
NO.data.station <- setdiff(ebc.ancillary$'station', union(names(ebc.icos.data), names(ebc.neon.data)))
if (length(NO.data.station) != 0) {
  ebc.ancillary <- ebc.ancillary[-which(ebc.ancillary$'station' %in% NO.data.station), ]
}
sites <- ebc.ancillary$'station'

# Set elaborated data structures

# . EBC input scenarios data series: 
EBC.base.data <- vector(mode = 'list', length = length(sites)); 
names(EBC.base.data) <- sites

# . New data list: Sph, Sbio, NETRAD corrections
new.data <- vector(mode = 'list', length = length(sites)); 
names(new.data) <- sites

# . NETRAD correction summary list (statistics of correction in % at each  site)
NETRAD.correction.summary <- vector(mode = 'list', length = length(sites)); 
names(NETRAD.correction.summary) <- sites


# base data preparation loop .....................................
i.bd <- 1
for (i.bd in 1:length(sites)) {
  
  # upload the current site data 
  cur.site <- sites[i.bd]
  cat(paste0('\n+++++++++++++ Making base dataset for ', cur.site, ' +++++++++++++\n'))
  
  # EBC data
  base.data <- data.frame()
  if (ebc.ancillary$'network'[i.bd] == 'ICOS') { # ICOS
    base.data <- ebc.icos.data[[cur.site]]  
  } else {                                         # NEON
    base.data <- ebc.neon.data[[cur.site]]  
  }
  base.data$'TIMESTAMP_START' <- as.character(base.data$'TIMESTAMP_START')
  base.data$'TIMESTAMP_END' <- as.character(base.data$'TIMESTAMP_END')
  
  # Fluxnet data (to merge)
  flxnt.tmp.df <- data.frame()
  if (ebc.ancillary$'network'[i.bd] == 'ICOS') { # ICOS
    flxnt.tmp.df <- data.frame(TIMESTAMP_START = as.character(flxnt.icos.data[[cur.site]]$'TIMESTAMP_START'),  # possibly add or remove fluxnet variables from here
                               GPP = flxnt.icos.data[[cur.site]][, grep('SE', grep('GPP', names(flxnt.icos.data[[cur.site]]), value = T), invert = T, value = T)],
                               SW_IN_F = flxnt.icos.data[[cur.site]]$'SW_IN_F',
                               TA_F = flxnt.icos.data[[cur.site]]$'TA_F',      
                               VPD_F = flxnt.icos.data[[cur.site]]$'VPD_F')    
  } else {                                         # NEON
    flxnt.tmp.df <- data.frame(TIMESTAMP_START = as.character(flxnt.neon.data[[cur.site]]$'TIMESTAMP_START'),  # possibly add or remove fluxnet variables from here
                               GPP = flxnt.neon.data[[cur.site]][, grep('SE', grep('GPP', names(flxnt.neon.data[[cur.site]]), value = T), invert = T, value = T)],
                               SW_IN_F = flxnt.neon.data[[cur.site]]$'SW_IN_F',
                               TA_F = flxnt.neon.data[[cur.site]]$'TA_F',      
                               VPD_F = flxnt.neon.data[[cur.site]]$'VPD_F')    
  }
  
  # Merging
  base.data <- merge(base.data,
                     flxnt.tmp.df,
                     by = 'TIMESTAMP_START', 
                     all.x = TRUE)
  
  rm(flxnt.tmp.df)
  
  if (save.newdata) {
    new.data[[cur.site]] <- data.frame(site = cur.site,
                                       TIMESTAMP_START = base.data$'TIMESTAMP_START')
  }
  
  # Photosynthesis energy consumption  -----------------------------------------------
  # Heat storage change associated with photosynthesis/respiration 
  # - [McGloin 2018; Moderow 2009]
  # J = −μ * GPP
  # μ is the specific energy for conversion due to photosynthesis (10.88*10^6 J per kg of CO2)
  # GPP is the gross primary productivity (kg CO2 m-2 s-1 where μmol of CO2 were converted to kg using the molar mass of CO2 = 0.044 kg of CO2 per mol)
  # Moderow 2009: JC adds energy to the storage and should be therefore positive during day and is zero during
  
  # GPP based 
  Spho <- c()
  if (!all(is.na(base.data$'GPP')) & !all(base.data$'GPP' == 0)) {
    Spho <- -base.data$'GPP' * umol2mol * CO2.mm * pecf
  } else {
    Spho <- rep(NA, nrow(base.data))
  }
  
  # Check day/night values according to Moderow 2009
  # summary(Spho) #[W m-2] # REF: Sp amounted to less than 10 W m−2
  
  if (save.newdata) { new.data[[cur.site]]$'Spho' <- Spho }
  
  # Biomass heat storage -----------------------------------------------------
  # Heat storage change in biomass [McGloin 2018; Kuticoff 2019; Moderow 2009] 
  # Jveg = mveg * cveg * ΔTveg / ΔT
  # mveg is the wet biomass (kg m−2): approximated by multiplying the dry biomass by a factor of 2
  # cveg is the canopy-specific heat capacity = 2958 J K−1 kg−1 
  # ΔTveg is the change in average biomass temperature (K): approximated by layer-weighted average air temperatures (calculated using air temperature profile measurements between the surface and top of the canopy) 
  # NOTE: the biomass is already converted in wet biomass during the Ancillary data production  
  
  # Average biomass from field campaigns (kg m−2) 
  bio.data.dyn <- data.frame()
  if (ebc.ancillary$'network'[ebc.ancillary$'station' == cur.site] == 'ICOS') {  
    bio.data.dyn <- bio.data.dyn.icos[[cur.site]]
  } else {
    bio.data.dyn <- bio.data.dyn.neon[[cur.site]]
  }
  
  # Interpolate the dynamic biomass data for each half-hour
  bio.wet.tmp <- data.frame()
  bio.wet.tmp <- data.frame(BIOMASS_DATE = base.data$'TIMESTAMP_START')
  
  # Drop possible duplicated values
  if (any(duplicated(bio.data.dyn$'BIOMASS_DATE'))) {bio.data.dyn <- bio.data.dyn[-which(duplicated(bio.data.dyn$'BIOMASS_DATE')), ]}
  
  # Add HHMM (arbitrarly)
  bio.data.dyn$'BIOMASS_DATE' <- paste0(bio.data.dyn$'BIOMASS_DATE', '1200')
  
  # If the biomass was measured during EBC dataset ok, otherwise use the median as a constant value
  if (max(as.numeric(bio.data.dyn$'BIOMASS_DATE')) > min(as.numeric(bio.wet.tmp$'BIOMASS_DATE'))) {
    bio.wet.tmp <- merge(
      bio.wet.tmp,
      bio.data.dyn,
      by = 'BIOMASS_DATE',
      all.x = TRUE
    )
    bio.wet.tmp$BIOMASS <- zoo::na.approx(bio.wet.tmp$BIOMASS, na.rm=F) # inner linear intepolation
    bio.wet.tmp$BIOMASS <- data.table::nafill(bio.wet.tmp$BIOMASS, type = 'nocb') # left fringe extrapolation
    bio.wet.tmp$BIOMASS <- data.table::nafill(bio.wet.tmp$BIOMASS, type = 'locf') # right fringe extrapolation
    
    bio.wet <- bio.wet.tmp$'BIOMASS'
    
    rm(bio.data.dyn, bio.wet.tmp)
    
  } else {
    
    bio.wet <- rep(ebc.ancillary$'bio_avg'[ebc.ancillary$'station' == cur.site], length(base.data$'TA_F')) 
    
  }
  
  
  # Average biomass temperature (K) ..........................
  # approximated by layer-weighted average air temperatures 
  
  TA.bio <- c()
  
  if (any(grepl(cur.site, TA.lw.files))) {
    TA.lw.tmp <- data.frame()
    TA.lw.tmp <- fread(grep(cur.site, TA.lw.files, value = TRUE), header = TRUE, sep=',', data.table = F)
    names(TA.lw.tmp)[1] <- 'TIMESTAMP_START'
    TA.lw.tmp$'TIMESTAMP_START' <- as.character(TA.lw.tmp$'TIMESTAMP_START')
    TA.lw.tmp$'TIMESTAMP_START' <- substr(TA.lw.tmp$'TIMESTAMP_START', 1, 12) # drop the seconds
    # Merge and align
    TA.lw.tmp <- merge(TA.lw.tmp,
                       data.frame(TIMESTAMP_START = base.data$'TIMESTAMP_START'),
                       by = 'TIMESTAMP_START',
                       all.y = TRUE)
    if (length(base.data$'TIMESTAMP_START') != length(TA.lw.tmp$'TA')) {
      stop('\nThe TA layer-weighted is not as long as the base.data. Check it out!\n')
    }
    
    TA.bio <- TA.lw.tmp$'TA'
    
    # fill the gaps with the gapfilled TA 
    if (any(is.na(TA.bio))) {
      TA.bio[is.na(TA.bio)] <- base.data$'TA_F'[is.na(TA.bio)]
    }
    
    rm(TA.lw.tmp)
    
  } else {
    
    # Se non è disponibile la TA layer-weighted, usa la TA @EC level
    if (!any(is.na(base.data$'TA_F'))) { # Take the gapfilled TA from the FLuxnet product
      TA.bio <- base.data$'TA_F'
    } else {
      if (all(is.na(base.data$'TA_F'))) { # fo "US-xGU" the fluxnet variables are all 0 or NA, take the TA from the fluxes data
        TA.bio <- zoo::na.approx(base.data$'TA', na.rm = FALSE)
      } else {
        TA.bio <- zoo::na.approx(base.data$'TA_F', na.rm = FALSE) 
      }
    }
    
  }
  
  # rough despiking
  if (any(na.omit(TA.bio) < -20 | na.omit(TA.bio) > 45)) {
    TA.bio[TA.bio < -20 | TA.bio > 45] <- NA
  }
  
  # convert to K
  if (median(TA.bio, na.rm = T) < 100) {
    TA.bio <- TA.bio + 273.14
  } 
  
  # change in average biomass temperature (K)
  dTbio <- c()
  dTbio <- c(NA, diff(TA.bio)) # add NA to replace the first record discarded by diff
  
  # rough despiking
  if (any(na.omit(dTbio) < -5 | na.omit(dTbio) > 5)) {
    dTbio[dTbio < -5 | dTbio > 5] <- NA 
  }
  
  # linear gapfilling (done now rather than later on Sbio, hopefully better)
  dTbio <- zoo::na.approx(dTbio, na.rm = FALSE)
  if (any(is.na(dTbio))) {
    dTbio <- zoo::na.locf(dTbio, na.rm = FALSE)
    dTbio <- zoo::na.locf0(dTbio, fromLast = TRUE)
  }
  
  # Heat storage change in biomass
  Sbio <- c()
  Sbio <- bio.wet * cshc * dTbio / 1800
  
  # Check it out
  # REF: Jveg ± 20 W m−2 [Moderow 2009]; Sc ± 5 W m−2 [Georg 2016]
  # summary(Sbio) #[W m-2]  
  
  # Rough despiking: drop the outliers out from the 0.01% and 99.99% quantiles 
  Qsbio <- quantile(Sbio, probs = c(0.01,0.99), na.rm = TRUE)
  Sbio[Sbio < Qsbio[1] | Sbio > Qsbio[2]] <- NA
  
  # linear gapfilling (done now rather than later on Sbio, hopefully better)
  Sbio <- zoo::na.approx(Sbio, na.rm = FALSE)
  if (any(is.na(Sbio))) {
    Sbio <- zoo::na.locf(Sbio, na.rm = FALSE)
    Sbio <- zoo::na.locf0(Sbio, fromLast = TRUE)
  }
  # summary(Sbio) #[W m-2]  
  
  if (save.newdata) {new.data[[cur.site]]$'Sbio' <- Sbio }
  
  rm(bio.wet, TA.bio, dTbio, Qsbio)
  
  
  
  # Net Radiation corrections --------------------------------------------------
  # NOTE_20250819: These variants of the NETRAD estimation must NOT be included in the AE improvemnt analysis,
  #                but ONLY to compare the EBC stats with those obtained using the original NETRAD.
  #                The aim is demonstrate their relevance for an actual EBC closure improvement
  
  # Daytime-nighttime splitting ................................................
  # Solar zenith angle
  sza_rad <- c()
  sza_rad <- solar_zenith(latitude = ebc.ancillary$'lat'[ebc.ancillary$'station' == cur.site], 
                          J = as.numeric(format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M', tz='GMT'), '%j')),  # DOY
                          solar_time_decimal = as.numeric(format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M', tz='GMT'), "%H")) + 
                            as.numeric(format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M', tz='GMT'), "%M")) / 60, # HoD as decimal at .5 resolution
                          sza.out='radians')
  
  # Set daytime-nightime based on solar elevation angle
  # FALSE durante la notte (incluso crepuscolo se sotto -0.833° se method = 'apparent')
  daynight_df <- data.frame(NULL)
  daynight_df <- classify_daynight(theta_z = sza_rad, 
                                   datetime = strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M', tz='GMT'),
                                   method = "apparent")
  daytime <- c(NULL)
  daytime <- daynight_df$'is_day_appar'
  
  
  ## 1) Slope correction ----------------------------------
  # Correct the measurement of shortwave solar radiation (incoming shortwave radiation) when instruments are mounted on inclined 
  # terrain, based on the method described by Olmo et al. (1999) and general considerations found in the works of Leuning et al. (2012) and McGloin et al. (2018).
  # The main idea is to calculate the corrected incident radiation on the inclined surface from the global radiation measured on
  # the horizontal plane, using the solar angles (zenith and azimuth angles), the slope angle of the terrain, and the aspect (orientation):
  
  # Albedo preliminary QC
  # Set albedo[%] == 0 during nightime 
  surface_albedo <- base.data$'ALB'
  # surface_albedo[cos(sza_rad) < 0] <- 0
  surface_albedo[!daytime] <- 0
  # linear interpolation
  surface_albedo <- zoo::na.approx(surface_albedo, na.rm = FALSE)
  
  # Load the station's DEM, if available
  if (ebc.ancillary$'network'[i.bd] == 'ICOS') { # The DEM correction applies only to ICOS data (for which we have the DEM)
    cur.site.dem <- grep(cur.site, list.files('C:/WORK/R/PROJECTS/EBC-ICOS/EBC_DATA/EBC_DTM_data/', full.names = TRUE), value = TRUE)
  } else {
    cur.site.dem <- NULL
  }
  
  # Compute the slope-corrected SW_IN
  # NOTE_20250810: questa versione di rad_topo_correction ha anche l'opzione per calcolare e correggere l'eventuale ombreggiamento da versante in base a DEM
  
  rad_topo_correction.df <- data.frame()
  rad_topo_correction.df <- rad_topo_correction(datetime = strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M', tz='GMT'), # POSIXct vector in UTC
                                                Gh = base.data$'SW_IN_F', #  numeric vector, global horizontal irradiance (W/m2)
                                                lat = ebc.ancillary$'lat'[ebc.ancillary$'station' == cur.site], #  station latitude (decimal degrees, N+)
                                                lon = ebc.ancillary$'lon'[ebc.ancillary$'station' == cur.site], #  station longitude (decimal degrees, E+)
                                                tilt = ebc.ancillary$'slope'[ebc.ancillary$'station' == cur.site], #  plane tilt (degrees) >> SLOPE
                                                azimuth_plane = ebc.ancillary$'aspect'[ebc.ancillary$'station' == cur.site], #  plane azimuth (degrees, 0=N, clockwise)  >> ASPECT
                                                albedo = surface_albedo, #  ground albedo
                                                add_reflection = TRUE, #  logical
                                                dem_corr = TRUE, # the DEM horizon mask: whether to force G values to 0 when the sensor is shadowed 
                                                dem = cur.site.dem, # optional, filename or SpatRaster of DEM
                                                max_dist = 5000,  # maximum search distance for horizon (m)
                                                az_step = 1, #  azimuth step for horizon calculation (degrees)
                                                sensor_height = ebc.ancillary$'EC_height'[ebc.ancillary$'station' == cur.site], #  sensor height above ground (m)
                                                night_control = FALSE #  whether to force nigh time G values to 0. # TRUE: force G to 0 during night, FALSE: leave the SW_IN untouched
  ) 

  # despiking
  SW_IN_slope <- c(NULL)
  SW_IN_slope <- rad_topo_correction.df$'G_tilt_total'
  SW_IN_slope[SW_IN_slope < - 50 | SW_IN_slope > 2000] <- NA
    
  # NETRAD correction on inclined plane .........................
  
  NETRAD_slope <- base.data$'NETRAD' - base.data$'SW_IN_F'
  NETRAD_slope <- NETRAD_slope + SW_IN_slope
    
  # Save statistics
  NRs.df <- data.frame(NRo=base.data$'NETRAD', NRs=NETRAD_slope)
  NRs.df <- na.omit(NRs.df)
  NRs.df[NRs.df == 0] <- NA
  NETRAD_slope.stats. <- round((quantile((NRs.df$'NRo' - NRs.df$'NRs') / NRs.df$'NRo', seq(0,1,0.1), na.rm=TRUE) *100), 2)
  NETRAD.correction.summary[[cur.site]][['NETRAD_slope.stats']] <- NETRAD_slope.stats.
  rm(NRs.df, NETRAD_slope.stats.)
  
  if (save.RAD_slope_correction & !all(is.na(SW_IN_slope))) {
    
    ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')
    
    # SWIN
    png(paste0('working_figures/', cur.site, '_rc11_slope_correction_SWIN.png'))
    # windows()
    sw.dc <- tapply(base.data$'SW_IN', ts.hh, mean, na.rm=T)
    sw.dc.se <- tapply(base.data$'SW_IN', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
    sw_s.dc <- tapply(SW_IN_slope, ts.hh, mean, na.rm=T)
    sw_s.dc.se <- tapply(SW_IN_slope, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
    plot(sw.dc, type='l', col=2, ylim=range(c(sw.dc, sw_s.dc)))
    polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(sw.dc + sw.dc.se), rev(sw.dc - sw.dc.se)),
            col=adjustcolor(8,0.5), border = NA)
    abline(v=24, lty=32)
    points(sw_s.dc, type='l', col=1)
    polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(sw_s.dc + sw_s.dc.se), rev(sw_s.dc - sw_s.dc.se)),
            col=adjustcolor(8,0.5), border = NA)
    legend('topleft', legend = c('SW_IN', 'SW_IN_slope'), lty=1, col=c(2,1), bty = 'n')
    mtext(cur.site, 3, 0)
    dev.off()
    
    # NETRAD
    png(paste0('working_figures/', cur.site, '_rc12_slope_correction_NETRAD.png'))
    # windows()
    netrad.dc <- tapply(base.data$'NETRAD', ts.hh, mean, na.rm=T)
    netrad.dc.se <- tapply(base.data$'NETRAD', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
    plot(netrad.dc, type='l', col=2)
    polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad.dc + netrad.dc.se), rev(netrad.dc - netrad.dc.se)),
            col=adjustcolor(2,0.5), border = NA)
    abline(v=24, lty=32)
    netrad_c.dc <- tapply(NETRAD_slope, ts.hh, mean, na.rm=T)
    netrad_c.dc.se <- tapply(NETRAD_slope, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
    points(netrad_c.dc, type='l', col=1)
    polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad_c.dc + netrad_c.dc.se), rev(netrad_c.dc - netrad_c.dc.se)),
            col=adjustcolor(8,0.5), border = NA)
    legend('topleft', legend = c('NETRAD', 'NETRAD_slope'), lty=1, col=c(2,1), bty = 'n')
    mtext(cur.site, 3, 0)
    dev.off()
  }
  
  if (save.newdata) { 
    new.data[[cur.site]]$'SW_IN_slope' <- SW_IN_slope
    new.data[[cur.site]]$'NETRAD_slope' <- NETRAD_slope 
  }
  
  
  ## 2) DEM shadowing correction ---------------------------------
  
  if (length(cur.site.dem) != 0) { # IF the DEM is available
    
    # Correzione SW_IN
    SW_IN_toposhadow <- c()
    
	# New method: SW_IN replaced by SW_DIFF when shadowed 
    
    # Diffuse radiation, estimated from SWIN if not measured
    SW_DIF <- c()
    if (any(grepl('SW_DIF', names(base.data)))) {
      
      SW_DIF <- base.data$'SW_DIF'
      
    } else {
      
      # Clearness index
      k_t <- rad_topo_correction.df$'kt'
      
      # condizioni validità 
      k_t.valid <- !is.na(k_t) &
        !is.na(base.data$'SW_IN_F') &
        # !is.na(SW_TOA) &
        # SW_TOA > 10 &     # evita crepuscolo instabile
        base.data$'SW_IN_F' > 5 # soglia x W m⁻²standard pratico SWIN f_diffuse = NA
      
      # Parametrizzazione Diffuse fraction (Reindl et al., 1990 / Erbs-like)
      
      f_diffuse <- rep(NA_real_, length(k_t))
      
      f_diffuse[k_t.valid & k_t <= 0.22] <- 1.0
      
      idx <- k_t.valid & k_t > 0.22 & k_t <= 0.8
      f_diffuse[idx] <- 1 - 1.09 * k_t[idx]
      
      f_diffuse[k_t.valid & k_t > 0.8] <- 0.18
      
      # sicurezza
      f_diffuse <- pmin(pmax(f_diffuse, 0), 1)
      
      # Stima radiazione diffusa
      SW_DIF <- base.data$'SW_IN_F' * f_diffuse
    }
    
    # Correzione
    is_shadowed <- !rad_topo_correction.df$'in_sun'

    SW_IN_toposhadow <- base.data$'SW_IN_F' # Inizializza con il valore non corretto
    SW_IN_toposhadow[is_shadowed & !is.na(SW_DIF)] <- SW_DIF[is_shadowed & !is.na(SW_DIF)]
    

    # NETRAD shadowing correction .........................
    NETRAD_toposhadow <- c()
    NETRAD_toposhadow <- base.data$'NETRAD' - base.data$'SW_IN_F'
    NETRAD_toposhadow <- NETRAD_toposhadow + SW_IN_toposhadow
    
    # Save statistics
    NRd.df <- data.frame(NRo=base.data$'NETRAD', NRd=NETRAD_toposhadow)
    NRd.df <- na.omit(NRd.df)
    NRd.df[NRd.df == 0] <- NA
    NETRAD_cur.site.demstats. <- round((quantile((NRd.df$'NRo' - NRd.df$'NRd') / NRd.df$'NRo', seq(0,1,0.1), na.rm=TRUE) *100), 2)
    NETRAD.correction.summary[[cur.site]][['NETRAD_toposhadow.stats']] <- NETRAD_cur.site.demstats.
    rm(NRd.df, NETRAD_cur.site.demstats.)
        
    if (save.RAD_dem_correction & !all(is.na(NETRAD_toposhadow))) {
      if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
      png(paste0('working_figures/', cur.site, '_rc20_toposhadow_correction_NETRAD.png'))
      # windows()
      netrad.dc <- tapply(base.data$'NETRAD', ts.hh, mean, na.rm=T)
      netrad.dc.se <- tapply(base.data$'NETRAD', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      plot(netrad.dc, type='l', col=2)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad.dc + netrad.dc.se), rev(netrad.dc - netrad.dc.se)),
              col=adjustcolor(2,0.5), border = NA)
      abline(v=24, lty=32)
      netrad_c.dc <- tapply(NETRAD_toposhadow, ts.hh, mean, na.rm=T)
      netrad_c.dc.se <- tapply(NETRAD_toposhadow, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      points(netrad_c.dc, type='l', col=1)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad_c.dc + netrad_c.dc.se), rev(netrad_c.dc - netrad_c.dc.se)),
              col=adjustcolor(8,0.5), border = NA)
      legend('topleft', legend = c('NETRAD', 'NETRAD_shadowing'), lty=1, col=c(2,1), bty = 'n')
      mtext(cur.site, 3, 0)
      dev.off()
      
    }
    
  } else { # IF the DEM is not available
    
    SW_IN_toposhadow <- rep(NA, length(base.data$'SW_IN_F'))
    NETRAD_toposhadow <- rep(NA, length(base.data$'NETRAD'))
    
  }
  
  
  if (save.newdata) { 
    new.data[[cur.site]]$'SW_IN_toposhadow' <- SW_IN_toposhadow
    new.data[[cur.site]]$'NETRAD_toposhadow' <- NETRAD_toposhadow 
  }
  
  
  ## 3) Tower disturbance correction ----
  # SW_OUT and LW_OUT disturbance due to the portion of the tower/pole falling within the radiometer's Field of View (FoV).
  # Radiometer mounted on the same truss (on an overhanging arm).
  # The "disturbance" is the contamination of the outgoing radiance (reflected/upwelled) entering the FoV because part of the truss section
  # falls within the cone of view.
  
  
  # check on radiometer setup input data
  rad.setup.ok <- is.na(c(rad.setup$'h_radiometro'[rad.setup$'Site code' == cur.site], 
                          rad.setup$'w_traliccio'[rad.setup$'Site code' == cur.site], 
                          rad.setup$'d_braccio'[rad.setup$'Site code' == cur.site]))
  
  if (ebc.ancillary$'network'[i.bd] == 'ICOS' & !any(rad.setup.ok)) {
    
    cat(paste('\n Tower disturbance correction for ', cur.site, '...'))
    
    # canopy heights (dynamic or static)
    cur.hc.df <- data.frame()
    if (!is.null(hc.data.dyn[[cur.site]]) && median(hc.data.dyn[[cur.site]]$'hc', na.rm = T) != 0) {
      # Import the dynamic canopy heights
      cur.hc.df <- data.frame(TIMESTAMP = base.data$'TIMESTAMP_START')
      cur.hc.df <- merge(cur.hc.df,
                         hc.data.dyn[[cur.site]],
                         by = 'TIMESTAMP',
                         # all.x=TRUE
                         all=TRUE            # use all=T to use canopy heights measured before the period covered by the base data (e.g. FR-FBn)
      )
      # Fill the possible fringe NAs with LOCB|NOCF
      if (any(is.na(cur.hc.df$'hc'))) {
        cur.hc.df$'hc' <- data.table::nafill(cur.hc.df$'hc', type = 'nocb') # left fringe extrapolation
        cur.hc.df$'hc' <- data.table::nafill(cur.hc.df$'hc', type = 'locf') # right fringe extrapolation
      }
      # Trim the hc data back to the period covered by the base data
      if (nrow(cur.hc.df) > nrow(base.data)) {
        t0hc <- which(cur.hc.df$'TIMESTAMP' == base.data$'TIMESTAMP_START'[1])
        tzhc <- which(cur.hc.df$'TIMESTAMP' == base.data$'TIMESTAMP_START'[nrow(base.data)])
        cur.hc.df <- cur.hc.df[t0hc:tzhc, ]
      }
      
    } else {
      # Import the static (average) canopy height
      cur.hc.df$'hc' <- ebc.ancillary$'hc_avg'[ebc.ancillary$'station' == cur.site]
    }
    
    
    ### SW_OUT correction ----
    # Kidston et la, 2010, Eqs 5:  Afov = (1-dfr)*Ag + dfr*At
    # A=albedo, subscripts g and t refer to the ground and tower respectively, 
    # dfr = disturbance fraction (% of the FoV occupied by the tower)
    
    # Compute the intersection among the radiometer SW_OUT FOV and the tower

    rad_tower_disturbance.df <- data.frame()
    rad_tower_disturbance.df <- disturbo_radiometro(
      fov_deg = rad.FoV_SWOUT,                                                     # Field of view del PYRANOMETRO  (gradi)
      h_radiometro = rad.setup$'h_radiometro'[rad.setup$'Site code' == cur.site],  # Altezza sensore dal suolo (m)
      w_traliccio = rad.setup$'w_traliccio_net'[rad.setup$'Site code' == cur.site],# Larghezza traliccio (m)
      d_braccio = rad.setup$'d_braccio'[rad.setup$'Site code' == cur.site],        # Distanza orizzontale radiometro-traliccio (m)
      h_c = cur.hc.df$'hc',                                                        # Vettore altezza copertura vegetale (m)
      albedo = albedo_tower,                                                       # Albedo del traliccio (0-1)
      irradianza = base.data$'SW_IN_F'                                             # Vettore irradianza incidente (W/m²)
    )
    
    # from tibble to dataframe
    rad_tower_disturbance.df <- as.data.frame(rad_tower_disturbance.df)
    
    if (save.RAD_towdist_omega) {
      png(paste0('working_figures/', cur.site, '_rc31_towdist_correction_SWOUT_geometry.png'), width = 480*4, height = 480*2, res = 130)
      # windows()
      sel.hh <- sample(which(daytime), 1) # random half-hour
      plot_disturbo_all(fov_deg = rad.FoV_SWOUT, 
                        h_radiometro = rad.setup$'h_radiometro'[rad.setup$'Site code' == cur.site], 
                        w_traliccio = rad.setup$'w_traliccio_net'[rad.setup$'Site code' == cur.site], 
                        d_braccio = rad.setup$'d_braccio'[rad.setup$'Site code' == cur.site], 
                        h_c = cur.hc.df$'hc'[sel.hh])
      
      dev.off()
    }
    
    # CORRECTION of SW_OUT: ....................................................
    
    # a. Disturbance fraction: fraction of the radiomenter FoV occupied by the tower (%)
    dfr.sw <- c(NULL)
    dfr.sw <- rad_tower_disturbance.df$'frazione_disturbo'

    # b. Measured albedo: actually measured albedo, affected by the tower
    albedo_fov <- surface_albedo/100 # 0-1
    
    # c. Albedo corrected (only "vegetation"): theoretical albedo, without the tower influence
    # Eq 5:  Afov = (1-dfr) * Ag + dfr * At  
    albedo_veg <- pmax(0, (albedo_fov - dfr.sw * albedo_tower) / (1 - dfr.sw))
    # NOTE: use pmax to avoid negative values when dfr.sw * albedo_tower > Afov
    
    # d. Correction coefficient [%]
    Acoef <- albedo_veg / albedo_fov
    # set to 1 (no correction) values 0 albedo (during night)
    Acoef[Acoef == 0] <- 1 
    Acoef[is.na(Acoef)] <- 1 
        
    # Correct SWOUT
    SW_OUT_towdist <- c(NULL)
    SW_OUT_towdist <- base.data$'SW_OUT' * Acoef
    
    # plot the results
    if (save.RAD_towdist_correction) {
      
      png(paste0('working_figures/', cur.site, '_rc32_towdist_correction_SWOUT.png'), width = 480*2, height = 480)
      # windows()
      par(mfrow=c(1,2), pty='s')
      
      hist(replace(Acoef, Acoef==1, NA), breaks=pretty(Acoef, 50), xlim=c(0,2),
           main='SWOUT correction factor')
      abline(v=1, col=1, lty=32)
      mtext(cur.site, 3)
      
      if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
      
      swout.dc <- tapply(base.data$'SW_OUT', ts.hh, mean, na.rm=T)
      swout.dc.se <- tapply(base.data$'SW_OUT', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      swout_c.dc <- tapply(SW_OUT_towdist, ts.hh, mean, na.rm=T)
      swout_c.dc.se <- tapply(SW_OUT_towdist, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      plot(swout.dc, type='l', col=1, ylim=range(c(swout.dc, swout_c.dc)))
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(swout.dc + swout.dc.se), rev(swout.dc - swout.dc.se)),
              col=adjustcolor(8,0.5), border = NA)
      abline(v=24, lty=32)
      points(swout_c.dc, type='l', col=2)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(swout_c.dc + swout_c.dc.se), rev(swout_c.dc - swout_c.dc.se)),
              col=adjustcolor(2,0.5), border = NA)
      legend('topleft', legend = c('measured', 'corrected'), lty=1, lwd=2, col=c(1,2), bty = 'n')
      mtext(paste('SW_OUT :', cur.site), 3, 0)
      dev.off()
      
    }
    
    
    ## LW_OUT correction --------------------------------
    # Kidston et la, 2010, Eqs 6:  Efov = (1-dfr)*Eg + dfr*Et
    # E = emissivity, subscripts g and t refer to the ground and tower respectively, dfr = disturbance fraction (% of the tower portion into the FoV)
    
    # Compute the intersection among the radiometer LW_OUT FOV and the tower
    rad_tower_disturbance.df <- data.frame()
    rad_tower_disturbance.df <- disturbo_radiometro(
      fov_deg = rad.FoV_LWOUT,                                                     # Field of view del radiometro (gradi)
      h_radiometro = rad.setup$'h_radiometro'[rad.setup$'Site code' == cur.site],  # Altezza sensore dal suolo (m)
      w_traliccio = rad.setup$'w_traliccio_net'[rad.setup$'Site code' == cur.site],# Larghezza traliccio (m)
      d_braccio = rad.setup$'d_braccio'[rad.setup$'Site code' == cur.site],        # Distanza orizzontale radiometro-traliccio (m)
      h_c = cur.hc.df$'hc',                                                        # Vettore altezza copertura vegetale (m)
      albedo = albedo_tower,                                                       # Albedo del traliccio (0-1)
      irradianza = base.data$'SW_IN_F'                                             # Vettore irradianza incidente (W/m²)
    )
    
    # from tibble to dataframe
    rad_tower_disturbance.df <- as.data.frame(rad_tower_disturbance.df)
    
    if (save.RAD_towdist_omega) {
      png(paste0('working_figures/', cur.site, '_rc33_towdist_correction_LWOUT_geometry.png'), width = 480*4, height = 480*2, res = 130)
      # windows()
      # sel.hh <- sample(which(daytime), 1) # random half-hour
      plot_disturbo_all(fov_deg = rad.FoV_LWOUT, 
                        h_radiometro = rad.setup$'h_radiometro'[rad.setup$'Site code' == cur.site], 
                        w_traliccio = rad.setup$'w_traliccio_net'[rad.setup$'Site code' == cur.site], 
                        d_braccio = rad.setup$'d_braccio'[rad.setup$'Site code' == cur.site], 
                        h_c = cur.hc.df$'hc'[sel.hh])
      dev.off()
    }
    
    # Surface emissivity computation .................
    
    # possibly fill NAs in the base.data$'TA_F' data
    if (any(is.na(base.data$'TA_F'))) {
      T_air <- base.data$'TA_F'
      T_air <- zoo::na.approx(T_air, na.rm=F) # inner linear intepolation
      T_air <- data.table::nafill(T_air, type = 'nocb') # left fringe extrapolation
      T_air <- data.table::nafill(T_air, type = 'locf') # right fringe extrapolation
    } else {
      T_air <- base.data$'TA_F'+ 273.14
    } 
    
    # From °C to K
    if (mean(T_air, na.rm=T) < 100) {
      T_air <- T_air + 273.15
    }
    
    surf_emiss.df <- data.frame()
    surf_emiss.df <- surf_emissivity(L_up_meas = base.data$'LW_OUT',              # W/m2   (misura upwelling)
                                     T_surface = NULL,                            # K o NULL
                                     T_air = T_air,                               # K (usato se T_surface NULL)
                                     deltaT_surface = ifelse(daytime,             # ΔT, correzione da applicare a T_air se T_surface mancante
                                                             2,                   #  +2 K di notte
                                                             4),                  #  +4 K di giorno
                                     simplified = FALSE,                          # se TRUE (sensore vicino a surface) usa la funzione semplice ε = L↑ / σTs^4, se FALSE (sensore misura a distanza) 
                                     L_down = replace(base.data$'LW_IN', is.na(base.data$'LW_IN'), 0),           # W/m2   (downwelling incident sulla superficie), scalare o vettore
                                     tau = 1,                                     # trasmittanza atmosfera (0-1), scalar o vettore
                                     L_path = 0,                                  # W/m2   emissione lungo la linea di nirta (scalar o vect)
                                     sigma_L = abs(base.data$'LW_OUT' * 0.1),     # errore su L_up_meas (W/m2), scalar o vect, opzionale: considerate specs CNR4: Pyrgeometer uncertainty	< 10% daily total [https://www.campbellsci.com/cnr4]
                                     sigma_T = 0.1+273.15                         # errore su T (K), scalar o vect, opzionale:  considerate specs HC2S3: ±0.1°C @ 23°C [https://www.campbellsci.eu/hc2s3#row360_]
    )
    
    if (save.RAD_towdist_emissivity) {
      
      emissivity_df <- data.frame()
      emissivity_df <- data.frame(
        TIMESTAMP=base.data$'TIMESTAMP_START',
        LW_OUT = base.data$'LW_OUT',
        emissivity = surf_emiss.df$'emissivita',
        e_sigma = surf_emiss.df$emissivita_sigma
      )
      
      library(tidyverse)
      library(lubridate)
      
      # --- dati e daily cycle come prima ---
      df2 <- emissivity_df %>%
        mutate(
          TIMESTAMP = ymd_hm(TIMESTAMP),
          month = month(TIMESTAMP, label = TRUE, abbr = TRUE),
          time_of_day = hour(TIMESTAMP) + minute(TIMESTAMP)/60
        )
      
      # --- Daily cycle mensile ---
      df_daily <- df2 %>%
        group_by(month, time_of_day) %>%
        summarise(
          LW_OUT     = mean(LW_OUT, na.rm = TRUE),
          emissivity = mean(emissivity, na.rm = TRUE),
          e_sigma    = mean(e_sigma, na.rm = TRUE),
          .groups = "drop"
        ) %>%
        mutate(
          LW_low  = LW_OUT * 0.90,
          LW_high = LW_OUT * 1.10,
          e_low   = emissivity - e_sigma/100,
          e_high  = emissivity + e_sigma/100
        )
      
      # --- Fissiamo scale ---
      LW_min <- 300
      LW_max <- 480
      emis_min <- 0
      emis_max <- 1.5
      
      # scale_factor <- (LW_max - LW_min) / (emis_max - emis_min)
      scale_factor <- 500
      
      df_daily <- df_daily %>%
        mutate(
          emissivity_scaled = pmin(pmax(emissivity * scale_factor, LW_min), LW_max),
          e_low_scaled      = pmin(pmax(e_low * scale_factor, LW_min), LW_max),
          e_high_scaled     = pmin(pmax(e_high * scale_factor, LW_min), LW_max)
        )
      
      # --- Plot corretto ---
      p <- ggplot(df_daily, aes(x = time_of_day)) +
        
        # Ribbon LW_OUT
        geom_ribbon(aes(ymin = LW_low, ymax = LW_high),
                    fill = "blue", alpha = 0.15) +
        
        # Ribbon emissività scalata
        geom_ribbon(aes(ymin = e_low_scaled, ymax = e_high_scaled),
                    fill = "red", alpha = 0.15) +
        
        # Linee
        geom_line(aes(y = LW_OUT), color = "blue", size = 0.8) +
        geom_line(aes(y = emissivity_scaled), color = "red", size = 0.8) +
        
        facet_wrap(~ month, ncol = 3) +
        
        scale_y_continuous(
          name = "LW_OUT (W/m²)",
          sec.axis = sec_axis(
            ~ . / scale_factor,
            name = "Emissivity",
            breaks = seq(0, 1.5, 0.3),
            labels = seq(0, 1.5, 0.3)
          )
        ) +
        scale_x_continuous(breaks = seq(0, 24, 3)) +
        labs(
          x = "Hour of day",
          title = "Monthly Mean Daily Cycles of LW_OUT and Emissivity"
        ) +
        theme_bw() +
        theme(
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 12),
          plot.title = element_text(size = 15),
          axis.title.y.left  = element_text(color = "blue"),
          axis.title.y.right = element_text(color = "red")
        ) +
        coord_cartesian(ylim = c(LW_min, LW_max))
      
      
      png(paste0('working_figures/', cur.site, '_rc34_towdist_correction_LWOUT_emissivity.png'), width = 480*4, height = 480*2, res = 130)
      print(p)
      dev.off()
      
      rm(df2, df_daily, scale_factor, p)
    }
    
    # CORRECTION of LW_OUT: ....................................................
    
    # a. Disturbance fraction: fraction of the radiomenter FoV occupied by the tower (%)
    dfr.lw <- c(NULL)
    dfr.lw <- rad_tower_disturbance.df$'frazione_disturbo'
    
    # b. "Measured" emissivity: estimated emissivity, affected by the tower
    emissivity_fov <- surf_emiss.df$'emissivita'
    
    # c. Emissivity corrected (only "vegetation"): theoretical emissivity, without the tower influence
    # Eq 6:  Efov = (1-dfr) * Eg + dfr * Et  
    emissivity_veg <- (emissivity_fov - dfr.lw * emissivity_tower) / (1-dfr.lw)
    emissivity_veg <- pmin(pmax(emissivity_veg, 0), 1) # limit to 0-1
    
    # d. Correction coefficient [%]
    Ecoef <- emissivity_veg/emissivity_fov
    Ecoef[Ecoef == 0] <- 1  # set the few 0 values to NA 
    Ecoef[Ecoef == 'Inf'] <- 1 # set to 1 (no correction) values divided by 0 (0 albedo during night data)
        
    # Correct LWOUT
    LW_OUT_towdist <- c(NULL)
    LW_OUT_towdist <- base.data$'LW_OUT' * Ecoef
    
    # plot the results 
    if (save.RAD_towdist_correction) {
      
      png(paste0('working_figures/', cur.site, '_rc35_towdist_correction_LWOUT.png'), width = 480*2, height = 480)
      # windows()
      par(mfrow=c(1,2), pty='s')
      
      hist(replace(Ecoef, Ecoef==1, NA), breaks=pretty(Ecoef, 50), xlim=c(0,2),
           main='SWOUT correction factor')
      abline(v=1, col=1, lty=32)
      mtext(cur.site, 3)
      
      if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
      
      lwout.dc <- tapply(base.data$'LW_OUT', ts.hh, mean, na.rm=T)
      lwout.dc.se <- tapply(base.data$'LW_OUT', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      lwout_c.dc <- tapply(LW_OUT_towdist, ts.hh, mean, na.rm=T)
      lwout_c.dc.se <- tapply(LW_OUT_towdist, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      plot(lwout.dc, type='l', col=1, ylim=range(c(lwout.dc, lwout_c.dc)))
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(lwout.dc + lwout.dc.se), rev(lwout.dc - lwout.dc.se)),
              col=adjustcolor(8,0.5), border = NA)
      abline(v=24, lty=32)
      points(lwout_c.dc, type='l', col=2)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(lwout_c.dc + lwout_c.dc.se), rev(lwout_c.dc - lwout_c.dc.se)),
              col=adjustcolor(2,0.5), border = NA)
      legend('topleft', legend = c('measured', 'corrected'), lty=1, lwd=2, col=c(1,2), bty = 'n')
      mtext(paste('LW_OUT :', cur.site), 3, 0)
      dev.off()
      
    }
    
    
    # NETRAD correction tower disturbance .........................
    
    NETRAD_towdist <- c()
    NETRAD_towdist <- base.data$'NETRAD' + base.data$'SW_OUT' + base.data$'LW_OUT'
    NETRAD_towdist <- NETRAD_towdist - SW_OUT_towdist - LW_OUT_towdist
    
    # Save statistics
    NRt.df <- data.frame(NRo=base.data$'NETRAD', NRt=NETRAD_towdist)
    NRt.df <- na.omit(NRt.df)
    NRt.df[NRt.df == 0] <- NA
    NETRAD_tower.stats. <- round((quantile((NRt.df$'NRo' - NRt.df$'NRt') / NRt.df$'NRo', seq(0,1,0.1), na.rm=TRUE) *100), 2)
    NETRAD.correction.summary[[cur.site]][['NETRAD_towdist.stats']] <- NETRAD_tower.stats.
    rm(NRt.df, NETRAD_tower.stats.)
        
    if (save.RAD_towdist_correction) {
      
      if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
      
      png(paste0('working_figures/', cur.site, '_rc36_towdist_correction_NETRAD.png'))
      # windows()
      netrad.dc <- tapply(base.data$'NETRAD', ts.hh, mean, na.rm=T)
      netrad.dc.se <- tapply(base.data$'NETRAD', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      netrad_c.dc <- tapply(NETRAD_towdist, ts.hh, mean, na.rm=T)
      netrad_c.dc.se <- tapply(NETRAD_towdist, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      plot(netrad.dc, type='l', col=1, ylim=range(c(netrad.dc, netrad_c.dc), na.rm = T))
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad.dc + netrad.dc.se), rev(netrad.dc - netrad.dc.se)),
              col=adjustcolor(8,0.5), border = NA)
      abline(v=24, lty=32)
      points(netrad_c.dc, type='l', col=2)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad_c.dc + netrad_c.dc.se), rev(netrad_c.dc - netrad_c.dc.se)),
              col=adjustcolor(2,0.5), border = NA)
      legend('topleft', legend = c('measured', 'corrected'), lty=1, lwd=2, col=c(1,2), bty = 'n')
      mtext(paste('NETRAD:', cur.site), 3, 0)
      dev.off()
    }
    
    cat('done.\n')
    
  } else { # NEON
    
    SW_OUT_towdist <- rep(NA, nrow(base.data))
    LW_OUT_towdist <- rep(NA, nrow(base.data))
    NETRAD_towdist <- rep(NA, nrow(base.data))
    
  }
  
  gc()
  
  if (save.newdata) { 
    new.data[[cur.site]]$'SW_OUT_towdist' <- SW_OUT_towdist
    new.data[[cur.site]]$'LW_OUT_towdist' <- LW_OUT_towdist 
    new.data[[cur.site]]$'NETRAD_towdist' <- NETRAD_towdist 
  }
  # summary(new.data[[cur.site]])
  
  
  ## 4) Radiometer FOV representativity correction ---------------------------------
  # use:
  # vis.ratio (VIS) = ratio based on visible bands, used to correct SW_OUT
  # nir.ratio (NIR) = ratio based on near infra-red bands, used to correct LW_OUT
  # or
  # tot.ratio (TOT) = joined ratio based on visible and near infra-red bands, used to correct both SW_OUT and LW_OUT
  
  if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
    
    cat(paste('\n FOV representativity correction for ', cur.site, '...'))
    
    # Import the data for the current site
    rad.ffp.fov.cur <- data.frame()
    rad.ffp.fov.cur <- rad.ffp.fov[rad.ffp.fov$'SITE_CODE' == cur.site, ]
    
    if (nrow(rad.ffp.fov.cur) > 0) {
      
      # build a dataframe with the time montly specification as YYYYmmddHHMM
      cur.fov.df <- data.frame()
      cur.fov.df <- data.frame(TIMESTAMP = paste0('2024', sprintf("%02d", 1:12), '151230'),
                               sw_ratio = rad.ffp.fov.cur$'vis_nir',
                               lw_ratio = rad.ffp.fov.cur$'surf_temp'
      )
      # despiking
      if (any(cur.fov.df$'sw_ratio' > 5, na.rm = T)) { cur.fov.df$'sw_ratio'[cur.fov.df$'sw_ratio' > 5] <- NA }
      if (any(cur.fov.df$'lw_ratio' > 5, na.rm = T)) { cur.fov.df$'lw_ratio'[cur.fov.df$'lw_ratio' > 5] <- NA }
      
      # Fill the possible NAs with linear|LOCB|NOCF (at monthly time scale)
      if (any(is.na(cur.fov.df$'sw_ratio'))) {
        cur.fov.df$'sw_ratio' <- zoo::na.approx(cur.fov.df$'sw_ratio', na.rm=F)           # inner linear intepolation
        cur.fov.df$'sw_ratio' <- data.table::nafill(cur.fov.df$'sw_ratio', type = 'nocb') # left fringe extrapolation
        cur.fov.df$'sw_ratio' <- data.table::nafill(cur.fov.df$'sw_ratio', type = 'locf') # right fringe extrapolation
      }
      if (any(is.na(cur.fov.df$'lw_ratio'))) {
        cur.fov.df$'lw_ratio' <- zoo::na.approx(cur.fov.df$'lw_ratio', na.rm=F)           # inner linear intepolation
        cur.fov.df$'lw_ratio' <- data.table::nafill(cur.fov.df$'lw_ratio', type = 'nocb') # left fringe extrapolation
        cur.fov.df$'lw_ratio' <- data.table::nafill(cur.fov.df$'lw_ratio', type = 'locf') # right fringe extrapolation
      }
      
      # Merge with the basedata
      cur.fov.df <- merge(data.frame(TIMESTAMP = base.data$'TIMESTAMP_START'),
                          cur.fov.df,
                          by = 'TIMESTAMP',
                          all.x=TRUE
      )
      # summary(cur.fov.df)
      
      # Fill the possible NAs with linear|LOCB|NOCF (at half-hourly time scale)
      if (any(is.na(cur.fov.df$'sw_ratio'))) {
        cur.fov.df$'sw_ratio' <- zoo::na.approx(cur.fov.df$'sw_ratio', na.rm=F)           # inner linear intepolation
        cur.fov.df$'sw_ratio' <- data.table::nafill(cur.fov.df$'sw_ratio', type = 'nocb') # left fringe extrapolation
        cur.fov.df$'sw_ratio' <- data.table::nafill(cur.fov.df$'sw_ratio', type = 'locf') # right fringe extrapolation
      }
      # summary(cur.fov.df$'sw_ratio')
      if (any(is.na(cur.fov.df$'lw_ratio'))) {
        cur.fov.df$'lw_ratio' <- zoo::na.approx(cur.fov.df$'lw_ratio', na.rm=F)           # inner linear intepolation
        cur.fov.df$'lw_ratio' <- data.table::nafill(cur.fov.df$'lw_ratio', type = 'nocb') # left fringe extrapolation
        cur.fov.df$'lw_ratio' <- data.table::nafill(cur.fov.df$'lw_ratio', type = 'locf') # right fringe extrapolation
      }
      
      # SW_OUT correction ......................................
      SW_OUT_fovffp <- c(NULL)
      SW_OUT_fovffp <- base.data$'SW_OUT' / cur.fov.df$'sw_ratio' # based on visible+NIR+SWIR bands ratio
      
      # plot the results
      if (save.RAD_fovffp_correction & !all(is.na(cur.fov.df$'sw_ratio'))) {
        png(paste0('working_figures/', cur.site, '_rc4_SWOUT_fovffp_correction.png'), width = 480*2, height = 480)
        # windows()
        par(mfrow=c(1,2), pty='s')
        hist(cur.fov.df$'sw_ratio', breaks=pretty(cur.fov.df$'sw_ratio', 10), xlim=c(0.5,1.5), main='SWOUT correction factor')
        abline(h=0, v=1)
        mtext(cur.site, 3)
        
        if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
        
        swout.dc <- tapply(base.data$'SW_OUT', ts.hh, mean, na.rm=T)
        swout.dc.se <- tapply(base.data$'SW_OUT', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
        swout_c.dc <- tapply(SW_OUT_fovffp, ts.hh, mean, na.rm=T)
        swout_c.dc.se <- tapply(SW_OUT_fovffp, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
        plot(swout.dc, type='l', col=2, ylim=range(c(swout.dc, swout_c.dc)))
        polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(swout.dc + swout.dc.se), rev(swout.dc - swout.dc.se)),
                col=adjustcolor(2,0.5), border = NA)
        abline(v=24, lty=32)
        points(swout_c.dc, type='l', col=1)
        polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(swout_c.dc + swout_c.dc.se), rev(swout_c.dc - swout_c.dc.se)),
                col=adjustcolor(8,0.5), border = NA)
        legend('topleft', legend = c('SW_OUT', 'SW_OUT_fovffp'), lty=1, col=c(2,1), bty = 'n')
        mtext(cur.site, 3, 0)
        dev.off()
      }
      
      
      # LW_OUT correction ......................................
      LW_OUT_fovffp <- c(NULL)
      LW_OUT_fovffp <- base.data$'LW_OUT' / cur.fov.df$'lw_ratio' # based on surface reflectance bands ratio

      # plot the results
      if (save.RAD_fovffp_correction& !all(is.na(cur.fov.df$'lw_ratio'))) {
        png(paste0('working_figures/', cur.site, '_rc4_LWOUT_fovffp_correction.png'), width = 480*2, height = 480)
        # windows()
        par(mfrow=c(1,2), pty='s')
        hist(cur.fov.df$'lw_ratio', breaks=pretty(cur.fov.df$'lw_ratio', 10), xlim=c(0.5,1.5), main='LWOUT correction factor')
        abline(h=0, v=1)
        mtext(cur.site, 3)
        
        if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
        
        lwout.dc <- tapply(base.data$'LW_OUT', ts.hh, mean, na.rm=T)
        lwout.dc.se <- tapply(base.data$'LW_OUT', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
        lwout_c.dc <- tapply(LW_OUT_fovffp, ts.hh, mean, na.rm=T)
        lwout_c.dc.se <- tapply(LW_OUT_fovffp, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
        plot(lwout.dc, type='l', col=2, ylim=range(c(lwout.dc, lwout_c.dc)))
        polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(lwout.dc + lwout.dc.se), rev(lwout.dc - lwout.dc.se)),
                col=adjustcolor(2,0.5), border = NA)
        abline(v=24, lty=32)
        points(lwout_c.dc, type='l', col=1)
        polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(lwout_c.dc + lwout_c.dc.se), rev(lwout_c.dc - lwout_c.dc.se)),
                col=adjustcolor(8,0.5), border = NA)
        legend('topleft', legend = c('LW_OUT', 'LW_OUT_fovffp'), lty=1, col=c(2,1), bty = 'n')
        mtext(cur.site, 3, 0)
        dev.off()
      }
      
      
      # NETRAD correction FOV representativity .........................
      
      NETRAD_fovffp <- base.data$'NETRAD' + base.data$'SW_OUT' + base.data$'LW_OUT'
      NETRAD_fovffp <- NETRAD_fovffp - SW_OUT_fovffp - LW_OUT_fovffp
      
      # Save statistics
      NRf.df <- data.frame(NRo=base.data$'NETRAD', NRf=NETRAD_fovffp)
      NRf.df <- na.omit(NRf.df)
      NRf.df[NRf.df == 0] <- NA
      NETRAD_fov.stats. <- round((quantile((NRf.df$'NRo' - NRf.df$'NRf') / NRf.df$'NRo', seq(0,1,0.1), na.rm=TRUE) *100), 2)
      NETRAD.correction.summary[[cur.site]][['NETRAD_fovffp.stats']] <- NETRAD_fov.stats.
      rm(NRf.df, NETRAD_fov.stats.)
      
      # # # check it out!
      if (save.RAD_fovffp_correction & !all(is.na(NETRAD_fovffp))) {
        if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
        png(paste0('working_figures/', cur.site, '_rc4_NETRAD_fovffp_correction_DC.png'))
        # windows()
        netrad.dc <- tapply(base.data$'NETRAD', ts.hh, mean, na.rm=T)
        netrad.dc.se <- tapply(base.data$'NETRAD', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
        netrad_c.dc <- tapply(NETRAD_fovffp, ts.hh, mean, na.rm=T)
        netrad_c.dc.se <- tapply(NETRAD_fovffp, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
        plot(netrad.dc, type='l', col=2, ylim=range(c(netrad.dc, netrad_c.dc)))
        polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad.dc + netrad.dc.se), rev(netrad.dc - netrad.dc.se)),
                col=adjustcolor(2,0.5), border = NA)
        abline(v=24, lty=32)
        points(netrad_c.dc, type='l', col=1)
        polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad_c.dc + netrad_c.dc.se), rev(netrad_c.dc - netrad_c.dc.se)),
                col=adjustcolor(8,0.5), border = NA)
        legend('topleft', legend = c('NETRAD', 'NETRAD_fovffp'), lty=1, col=c(2,1), bty = 'n')
        mtext(cur.site, 3, 0)
        dev.off()
      }
      
      cat('done.\n')
      
    } else { # no data
      cat(paste('\n NO FOV representativity correction data for ', cur.site))
      SW_OUT_fovffp <- rep(NA, nrow(base.data))
      LW_OUT_fovffp <- rep(NA, nrow(base.data))
      NETRAD_fovffp <- rep(NA, nrow(base.data))
    }
    
  } else { # NEON
    
    SW_OUT_fovffp <- rep(NA, nrow(base.data))
    LW_OUT_fovffp <- rep(NA, nrow(base.data))
    NETRAD_fovffp <- rep(NA, nrow(base.data))
    
  }
  
  if (save.newdata) { 
    new.data[[cur.site]]$'SW_OUT_fovffp' <- SW_OUT_fovffp
    new.data[[cur.site]]$'LW_OUT_fovffp' <- LW_OUT_fovffp 
    new.data[[cur.site]]$'NETRAD_fovffp' <- NETRAD_fovffp 
  }
  # summary(new.data[[cur.site]])
  
  
  ## 5) FULL correction ---------------------
  # NOTE: note used for the EBC paper
  
  if (ebc.ancillary$'network'[i.bd] == 'ICOS'  & !any(rad.setup.ok) ) {
    
    cat(paste('\n RAD full correction for ', cur.site, '...'))
    
    NETRAD_fullcorr <- base.data$'NETRAD'
    
    # Slope correction .........................................................
    # use the corrected data directly
    NETRAD_fullcorr <- NETRAD_fullcorr - base.data$'SW_IN_F'
    NETRAD_fullcorr <- NETRAD_fullcorr + SW_IN_slope
    
    # Topographic shadowing correction .........................................
    # RE-apply the correction to the slope corrected values
    SW_IN_toposhadow.2 <- c()
    SW_IN_toposhadow.2 <- replace(SW_IN_slope, rad_topo_correction.df$'in_sun' == FALSE, 0)
    NETRAD_fullcorr <- NETRAD_fullcorr - SW_IN_slope
    NETRAD_fullcorr <- NETRAD_fullcorr + SW_IN_toposhadow.2
    
    # Tower disturbance correction .............................................
    # Can apply directly the sums as this correction works on OUTGOING values
    NETRAD_fullcorr <- NETRAD_fullcorr + base.data$'SW_OUT' + base.data$'LW_OUT'
    NETRAD_fullcorr <- NETRAD_fullcorr - SW_OUT_towdist - LW_OUT_towdist
    
    # Radiometer FOV representativity correction ...............................
    # Apply the fov correction to the tower disturbance corrected SWOUT and LWOUT
    SW_OUT_fovffp.2 <- SW_OUT_towdist / cur.fov.df$'sw_ratio'
    LW_OUT_fovffp.2 <- LW_OUT_towdist / cur.fov.df$'lw_ratio'
    NETRAD_fullcorr <- NETRAD_fullcorr + SW_OUT_towdist + LW_OUT_towdist
    NETRAD_fullcorr <- NETRAD_fullcorr - SW_OUT_fovffp.2 - LW_OUT_fovffp.2
    
    
    # Save statistics
    NRf.df <- data.frame(NRo=base.data$'NETRAD', NRf=NETRAD_fullcorr)
    NRf.df <- na.omit(NRf.df)
    NRf.df[NRf.df == 0] <- NA
    NETRAD_full.stats. <- round((quantile((
      NRf.df$'NRo' - NRf.df$'NRf') / NRf.df$'NRo', 
      seq(0,1,0.1), na.rm=TRUE) *100), 2)
    NETRAD.correction.summary[[cur.site]][['NETRAD_fullcorr.stats']] <- NETRAD_full.stats.
    rm(NRf.df, NETRAD_full.stats.)
    
    # Save diurnal cycle
    if (save.NETRAD_full_correction) {
      if (!'ts.hh' %in% ls()) {ts.hh <- format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M'), '%H:%M')}
      png(paste0('working_figures/', cur.site, 'rc5_NETRAD_fullcorrections_DC.png'))
      # windows()
      netrad.dc <- tapply(base.data$'NETRAD', ts.hh, mean, na.rm=T)
      netrad.dc.se <- tapply(base.data$'NETRAD', ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      netrad_c.dc <- tapply(NETRAD_fullcorr, ts.hh, mean, na.rm=T)
      netrad_c.dc.se <- tapply(NETRAD_fullcorr, ts.hh, function(x) sd(x, na.rm = T) / sqrt(length(x)))
      plot(netrad.dc, type='l', col=2)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad.dc + netrad.dc.se), rev(netrad.dc - netrad.dc.se)),
              col=adjustcolor(2,0.5), border = NA)
      points(netrad.dc, type='l', col=2, lwd=2)
      abline(v=24, lty=32)
      polygon(c(seq_along(unique(ts.hh)), rev(seq_along(unique(ts.hh)))), c(c(netrad_c.dc + netrad_c.dc.se), rev(netrad_c.dc - netrad_c.dc.se)),
              col=adjustcolor(4,0.5), border = NA)
      points(netrad_c.dc, type='l', col=4, lwd=2)
      legend('topleft', legend = c('NETRAD', 'NETRAD_FullCorr'), lty=1, col=c(2,4), bty = 'n')
      mtext(cur.site, 3, 0)
      dev.off()
    }
    
    cat('done!\n')
    
  } else { # NEON
    
    NETRAD_fullcorr <- rep(NA, nrow(base.data))
    
  }
  
  if (save.newdata) { 
    new.data[[cur.site]]$'NETRAD_fullcorr' <- NETRAD_fullcorr 
  }
  
  
  
  # EBC INPUTs scenarios ----------------------------------
  
  # Select the proper one-point storage fluxes:
  # in ICOS the station without profile have SH/SLE == SH_1P/SLE_1P (SH_1P/SLE_1P columns are not present at all):
  # se hanno profile S1P è vero, se no SH == S1P
  # in NEON all station have a profile, but still the 1P is computed
  if(cur.site == 'FR-Lus'){base.data$'SH_FLAG' <- 1} # fr-lus è una stazione processata con il codice vecchio e non ha SH_FLAG
  sh1p <- sle1p <- c(NULL)
  if (any(grepl('SH_FLAG', names(base.data)))) { # ICOS
    sh1p <- switch(all(base.data$'SH_FLAG' != 0, na.rm = TRUE)+1, base.data$'SH_1P', base.data$'SH') # se no profile si chiama SH ma è SH_1P
    sle1p <- switch(all(base.data$'SLE_FLAG' != 0, na.rm = TRUE)+1, base.data$'SLE_1P', base.data$'SLE')
  } else {                                # NEON
    sh1p <- base.data$'SH_1P'
    sle1p <- base.data$'SH_1P'
  }
  
  
  ## Available energy --------------------
  
  AE.df <- data.frame(NULL)
  AE.df <- data.frame(
    # Net radiation
    NR = base.data$'NETRAD',
    # Net radiation slope corrected
    NRslope = NETRAD_slope,
    # Net radiation DEM topo-shadow corrected
    NRdem = NETRAD_toposhadow,
    # Net radiation tower disturbance corrected
    NRtow = NETRAD_towdist,
    # Net radiation Rad_FOV corrected
    NRfov = NETRAD_fovffp,
    # Net radiation fully corrected
    NRfullcor = NETRAD_fullcorr,
    # Soil heat flux
    G2 = -base.data$'G_2',
    # Soil heat flux *storage corrected*  (same as base.data$'G_1' +  base.data$'SG')
    G1 = -base.data$'G_1',
    # One-point sensible heat storage fluxes (SC_FLAG: 0=profile, 1=1P, 2=LinInt, 3=RFlux::cleansetQC > if only 1P storage is present SH == SH_1P in the base data)
    SH1P = -sh1p,
    # 1 point latent heat storage fluxes
    SLE1P = -sle1p,
    # Profile sensible heat storage fluxes >> if only 1P storage is present SH is SH_1P
    SH = -base.data$'SH',
    # profile latent heat storage fluxes >> if only 1P storage is present SLE is SLE_1P
    SLE = -base.data$'SLE',
    # Photosynthesis energy consumption heat storage
    Spho = switch((!all(is.na(Spho))) + 1, rep(NA, length(base.data$'NETRAD')), -Spho),
    # Biomass heat storage
    Sbio = switch((!all(is.na(Sbio))) + 1, rep(NA, length(base.data$'NETRAD')), -Sbio)
    
  )

  ## Turbulent energy --------------------
  
  # Prepare single variables according to the processing keys
  
  # . Sensible and Latent heat fluxes UNCLEANED (just a rough outlier removal)
  h.un <- base.data$'H_UNCLEANED'; # h.un[h.un < -500 | h.un > 1000] <- NA # NOTE_20250908: this filter has been moved to DoInputData.R
  le.un <- base.data$'LE_UNCLEANED'; # le.un[le.un < -500 | le.un > 1000] <- NA  # NOTE_20250908: this filter has been moved to DoInputData.R
  
  # . Sensible and Latent heat fluxes QC FKN04 [0,1] (discard H_SSITC_TEST = 2) - ONLY ICOS DATA
  if ('H_SSITC_TEST' %in% names(base.data)) {
    h.fk01 <- h.un; h.fk01[base.data$'H_SSITC_TEST' == 2] <- NA
    le.fk01 <- le.un; le.fk01[base.data$'LE_SSITC_TEST' == 2] <- NA
  } else {
    h.fk01 <- rep(NA, nrow(base.data))
    le.fk01 <- rep(NA, nrow(base.data))
  }
  
  # . Sensible and Latent heat fluxes QC FKN04 [0] (discard H_SSITC_TEST = 1|2) - ONLY ICOS DATA
  # NOTE_20250819: the EBC regression results for this TE NOT need to go into the processing improvement plot 
  #                can be compared with the regression using the f00 fluxes (above)
  if ('H_SSITC_TEST' %in% names(base.data)) {
    h.fk00 <- h.un; h.fk00[base.data$'H_SSITC_TEST' > 0] <- NA
    le.fk00 <- le.un; le.fk00[base.data$'LE_SSITC_TEST' > 0] <- NA
  } else {
    h.fk00 <- rep(NA, nrow(base.data))
    le.fk00 <- rep(NA, nrow(base.data))
  }
  
  # . Sensible and Latent heat fluxes full QC [QC ICOS Vitale || NEON approach] AND in-situ spectral assessment
  # added in the dataframe below directly
  
  # . Sensible and Latent heat fluxes ANLYTICAL spectral assessment - ONLY ICOS DATA
  # taken from NRT data (full 2024)
  if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
    
    if (nrow(nrt.icos.data[[cur.site]]) != 0) {
      
      Hspanal.df <- data.frame()
      Hspanal.df <- data.frame(TIMESTAMP_START = nrt.icos.data[[cur.site]]$'TIMESTAMP_START',
                               Hspanal = nrt.icos.data[[cur.site]]$'H')
      Hspanal.df <- merge(x = Hspanal.df, y = data.frame(TIMESTAMP_START = base.data$'TIMESTAMP_START'),
                          by = 'TIMESTAMP_START',
                          all.y = TRUE)
      h.spanal <- Hspanal.df$'Hspanal'
      
      LEspanal.df <- data.frame()
      LEspanal.df <- data.frame(TIMESTAMP_START = nrt.icos.data[[cur.site]]$'TIMESTAMP_START',
                                LEspanal = nrt.icos.data[[cur.site]]$'LE')
      LEspanal.df <- merge(x = LEspanal.df, y = data.frame(TIMESTAMP_START = base.data$'TIMESTAMP_START'),
                           by = 'TIMESTAMP_START',
                           all.y = TRUE)
      le.spanal <- LEspanal.df$'LEspanal'
      
      rm(Hspanal.df, LEspanal.df)
      
      
    } else {
      
      h.spanal <- rep(NA, nrow(base.data))
      le.spanal <- rep(NA, nrow(base.data))
      
    }
    
  } else {
    
    h.spanal <- rep(NA, nrow(base.data))
    le.spanal <- rep(NA, nrow(base.data))
    
  }
  
  # . Sensible and Latent heat fluxes full QC + u* filtering
  h.ust <- base.data$'H'
  le.ust <- base.data$'LE'
  if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
    h.ust[base.data$'USTAR' < flxnt.icos.usttrhd[[cur.site]]] <- NA
    le.ust[base.data$'USTAR' < flxnt.icos.usttrhd[[cur.site]]] <- NA
  } else {
    h.ust[base.data$'USTAR' < flxnt.neon.usttrhd[[cur.site]]] <- NA
    le.ust[base.data$'USTAR' < flxnt.neon.usttrhd[[cur.site]]] <- NA
  }
  
  # # . Sensible and Latent heat fluxes full QC + footprint filtering  - ONLY ICOS DATA
  # # # NOTE_20250819: considered no more as it has an additional trivial contribution 
  # h.ffp <- base.data$'H'
  # le.ffp <- base.data$'LE'
  # if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
  #   h.ffp[base.data$'FOOTPRINT_TA_CONTR' < ffp.thrhld] <- NA
  #   le.ffp[base.data$'FOOTPRINT_TA_CONTR' < ffp.thrhld] <- NA
  # } else {
  #   h.ffp <- rep(NA, nrow(base.data))
  #   le.ffp <- rep(NA, nrow(base.data))
  # }
  
  # # . Sensible and Latent heat fluxes full QC + u* filtering + footprint filtering  - ONLY ICOS DATA
  # # # NOTE_20250819: considered no more as it has an additional trivial contribution
  # h.ustffp <- h.ust
  # le.ustffp <- le.ust
  # if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
  #   h.ustffp[base.data$'FOOTPRINT_TA_CONTR' < ffp.thrhld] <- NA
  #   le.ustffp[base.data$'FOOTPRINT_TA_CONTR' < ffp.thrhld] <- NA
  # }
  
  # Fill the TE data frame
  TE.df <- data.frame(NULL)
  TE.df <- data.frame(
    # Sensible and Latent heat fluxes UNCLEANED (just a rough outlier removal)
    Hun = h.un,
    LEun = le.un,
    # Sensible and Latent heat fluxes QC FKN04 [0,1] (discard H_SSITC_TEST = 2) - ONLY ICOS DATA
    Hfk01 = h.fk01,
    LEfk01 = le.fk01,
    # Sensible and Latent heat fluxes QC FKN04 [0] (discard H_SSITC_TEST = 1|2) - ONLY ICOS DATA
    Hfk00 = h.fk00,
    LEfk00 = le.fk00,
    # . Sensible and Latent heat fluxes full QC [QC ICOS Vitale || NEON approach] AND in-situ spectral assessment
    # NOTE: For ICOS these data must be used as reference for the spectral correction of Fratini 2012
    H = base.data$'H',
    LE = base.data$'LE',
    # Sensible and Latent heat fluxes spanalYTICAL spectral assessment - ONLY ICOS DATA
    Hspanal = h.spanal,
    LEspanal = le.spanal,
    # Sensible and Latent heat fluxes full QC + u* filtering
    Hust = h.ust,
    LEust = le.ust
    # # . Sensible and Latent heat fluxes full QC + footprint filtering  - ONLY ICOS DATA
    # Hffp = h.ffp,
    # LEffp = le.ffp,
    # # . Sensible and Latent heat fluxes full QC + u* filtering + footprint filtering  - ONLY ICOS DATA
    # Hustffp = h.ustffp,
    # LEustffp = le.ustffp
  )
  
  # # Check spectral effect .....
  if (save.spcomp & !(all(is.na(h.spanal)))) {
    png(paste0('working_figures/', cur.site, '_SpectralEffectComparison.png'), width = 24, height = 16, units = 'cm', res = 100)
    # windows()
    par(mfrow=c(2,3), mar=c(2,2,2,2), oma=c(2,2,0.1,0.1))
    # H
    plot(base.data$'H', h.spanal)
    mtext('in-situ vs analytical', 3, 0)
    abline(0,1)
    Hdiff <- base.data$'H' - h.spanal
    plot(Hdiff)
    mtext('in-situ - analytical', 3, 0)
    # Boxplot della variabile per classi di RH
    RH_class <- cut(base.data$'RH', breaks=seq(0, 100, by=5), right=FALSE, include.lowest=TRUE)
    boxplot(Hdiff ~ RH_class, xlab="Classi di RH (%)", ylab="Hdiff", las=2, cex.axis=0.7, outline = FALSE)
    abline(h=0)
    mtext('in-situ - analytical by RH', 3, 0)
    # LE
    plot(base.data$'LE', le.spanal)
    abline(0,1)
    LEdiff <- base.data$'LE' - le.spanal
    plot(LEdiff)
    # Boxplot della variabile per classi di RH
    # RH_class <- cut(base.data$'RH', breaks=seq(0, 100, by=5), right=FALSE, include.lowest=TRUE)
    boxplot(LEdiff ~ RH_class, xlab="Classi di RH (%)", ylab="LEdiff", las=2, cex.axis=0.7, outline = FALSE)
    abline(h=0)
    dev.off()
  }
  
  # free some memory
  rm(h.un,le.un,h.fk01, le.fk01, h.fk00, le.fk00, h.spanal, le.spanal, h.ust, le.ust) #h.ffp, le.ffp, h.ustffp, le.ustffp)
  
  
  ## Available energy COMBO --------------------
  # dataframe con tutti gli scenari di aggregazioni possibili
  # sommo fra loro le colonne del dataframe AE.df considerando tutte le possibili combinazioni di colonne (dalla somma di una sola colonna fino alla somma di tutte insieme).
  # creando un nuovo dataframe con una colonna per ogni combinazione possibile, in cui ogni colonna è la somma delle colonne corrispondenti di AE.df
  # filtrarlo per casi impossibili, tipo
  # . combinazioni con più di una netrad
  # . combinazioni con più di un G
  # . combinazioni con più di uno storage per lo stesso scalare (es: SH_1P e SH)
  # UPDATE_20250929: now all of this is implemented into a function ("make_combos_filtered").
  #                  if other cases need to be added/removed, act on the function
  

  # t0 <- Sys.time()
  AEcombo_df <- data.frame()
  cat('\nComputing AE combos ...')
  AEcombo_df <- make_combos_filtered(AE.df)
  # tz <- Sys.time()
  # tz-t0
  cat('done.\n')
  
  ### Relevant cases fetching ----
  # Bring forward ONLY the main cases, which have sense in the real world (for sake of simplicity in the interpretation of results)
  # names(AEcombo_df)
  
  # level_0: NETRAD ONLY
  l0.combos <- names(AEcombo_df)[grepl("NR$|NRslope$|NRdem$|NRtow$|NRfov$|NRfullcor$", names(AEcombo_df))]
  # level_1: NETRAD* + G2
  l1.combos <- names(AEcombo_df)[apply(sapply(c("_G2$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  # level_2: NETRAD* + G1
  l2.combos <- names(AEcombo_df)[apply(sapply(c("_G1$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  # level_3: NETRAD* + G1 + Sx1P
  l3.combos <- names(AEcombo_df)[apply(sapply(c("_G1_SH1P_SLE1P$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  # level_4: NETRAD* + G1 + Sx
  l4.combos <- names(AEcombo_df)[apply(sapply(c("_G1_SH_SLE$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  # level_5: NETRAD* + G1 + Sx + Spho
  l5.combos <- names(AEcombo_df)[apply(sapply(c("_G1_SH_SLE_Spho$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  # level_6: NETRAD* + G1 + Sx + Spho + Sbio
  l6.combos <- names(AEcombo_df)[apply(sapply(c("_G1_SH_SLE_Spho_Sbio$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  # level_6p: NETRAD_OR + G1 + Sx1P + Spho + Sbio (used for only demonstrative purpose in the analysis)
  l6p.combos <- names(AEcombo_df)[apply(sapply(c("NR_G1_SH1P_SLE1P_Spho_Sbio$"), function(p) grepl(p, names(AEcombo_df))), 1, all)]
  
  AEcombo_df <- AEcombo_df[, c(l0.combos, l1.combos, l2.combos, l3.combos, l4.combos, l5.combos, l6.combos, l6p.combos)]
  # names(AEcombo_df)
  
  # cat(paste('\nAfter the relevant case fetching,', ncol(AEcombo_df), 'AE combo remain for the EBC analysis.\n'))
  
  # free some memory
  rm(l0.combos, l1.combos, l2.combos, l3.combos, l4.combos, l5.combos, l6.combos, l6p.combos)
  
  
  ## Turbulent energy COMBO --------------------
  
  cat('\nComputing TE combos ...')
  
  # Estrai i nomi delle colonne H e LE
  cols_H <- colnames(TE.df)[grepl("^H", colnames(TE.df))]
  cols_LE <- colnames(TE.df)[grepl("^LE", colnames(TE.df))]
  
  # Lista per salvare i risultati
  TEcombo_list <- list()
  
  # Ciclo su ogni colonna H e LE singolarmente
  for (h_col in cols_H) {
    for (le_col in cols_LE) {
      # Somma riga per riga la coppia di colonne
      combined_sum <- TE.df[[h_col]] + TE.df[[le_col]]
      # Crea un nome che indica le colonne sommate
      col_name <- paste(h_col, le_col, sep = "_")
      # Salva nella lista
      TEcombo_list[[col_name]] <- combined_sum
    }
  }
  # Combina in dataframe finale
  TEcombo_df <- as.data.frame(TEcombo_list)
  
  rm(TEcombo_list, cols_H, cols_LE, combined_sum, col_name)
  gc()
  
  
  ### Impossible cases filtering ----
  
  # > combinazioni con H e LE filtrati con approcci diversi
  # Pattern per prefissi "H" e "LE" seguito da suffissi alfanumerici
  h_match <- str_extract(names(TEcombo_df), "H[a-z0-9]+|H")
  le_match <- str_extract(names(TEcombo_df), "LE[a-z0-9]+|LE")
  # Verifica che entrambi siano presenti
  if (!all(is.na(h_match)) | !all(is.na(le_match))) {
    # Rimuove il prefisso per ottenere solo i suffissi
    h_suff <- sub("^H", "", h_match)
    le_suff <- sub("^LE", "", le_match)
    # nomi con suffissi uguali
    same_hle <- na.omit(names(TEcombo_df)[h_suff == le_suff])
    # filtro
    TEcombo_df <- TEcombo_df[, same_hle]
    
    rm(h_suff, le_suff, same_hle,h_match, le_match)
  }
  
  
  # cat(paste('\nAfter the case removal,', ncol(TEcombo_df), 'TE combo remain for the EBC analysis.\n'))
  
  
  ### Relevant cases fetching ----
  # Bring forward ONLY the main cases, which have sense in the real world (for sake of simplicity in the interpretation of results)
  # >> all relevant cases are there already
  cat('done.\n')
  
  # free some memory
  gc()
  
  
  # GAPFILLED DATA ------------------------
  # NOTE_20250820: To use gapfilled data, AEcombo and TE combo dataframes are done directly, without the full-factorial sum (as for the non-gapfilled data above)
  #                Then the EBC analysis is computed separately, also adding the daily sums
  
  if (cur.site %in% c(names(ebc.GF.icos.data), names(ebc.GF.neon.data))) {
    
    # set -9999 values in ebc.GF.icos.data to NA ...
    tmp.GF.df0 <- data.frame()
    if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
      tmp.GF.df0 <- ebc.GF.icos.data[[cur.site]]
    } else {
      tmp.GF.df0 <- ebc.GF.neon.data[[cur.site]]
    }
    
    # Elimino i -9999: si può fare perchè sono solo in NR, solo nella prima parte dell'anno (che comunque va buttata)
    tmp.GF.df0[tmp.GF.df0 == '-9999'] <- NA
    
    # Align gapfilled data with the base data
    GF.df <- merge(
      data.frame(TIMESTAMP_END = base.data$'TIMESTAMP_END'),
      tmp.GF.df0,
      by = 'TIMESTAMP_END',
      all.x = TRUE
    )
    rm(tmp.GF.df0)
    
    
    ## Available energy --------------------
    AEcombo_GF_df <- data.frame()
    
    # level_0: NETRAD ONLY
    # Rough filtering on ONEFlux GF FLAG == 3 
    # !! NOTA: si può fare perchè sono solo nella prima parte dell'anno 
    GF.df$'NETRAD_FILLED'[1:30000][GF.df$'NETRAD_GF_QC'[1:30000] == 3] <- NA
    
    AEcombo_GF_df <- data.frame(
      NR = GF.df$'NETRAD_FILLED'
    )
    
    # . level_1: NETRAD + G2
    # Rough filtering on ONEFlux GF FLAG == 3 
    GF.df$'G_2_FILLED'[1:30000][GF.df$'G_2_GF_QC'[1:30000] == 3] <- NA
    
    AEcombo_GF_df$'NR_G2' <- AEcombo_GF_df$'NR' - GF.df$'G_2_FILLED'
    
    # . level_2: NETRAD + G1
    # Rough filtering on ONEFlux GF FLAG == 3 
    GF.df$'G_1_FILLED'[1:30000][GF.df$'G_1_GF_QC'[1:30000] == 3] <- NA
    
    AEcombo_GF_df$'NR_G1' <- AEcombo_GF_df$'NR' - GF.df$'G_1_FILLED'
    
    # . level_3: NETRAD + G1 + Sx1P
    # Rough filtering on ONEFlux GF FLAG == 3 
    GF.df$'SH_1P_FILLED'[1:30000][GF.df$'SH_1P_GF_QC'[1:30000] == 3] <- NA
    GF.df$'SLE_1P_FILLED'[1:30000][GF.df$'SLE_1P_GF_QC'[1:30000] == 3] <- NA
    
    AEcombo_GF_df$'NR_G1_SH1P_SLE1P' <- AEcombo_GF_df$'NR_G1' - GF.df$'SH_1P_FILLED' - GF.df$'SLE_1P_FILLED'
    
    # . level_4: NETRAD + G1 + Sx
    # Rough filtering on ONEFlux GF FLAG == 3 
    GF.df$'SH_FILLED'[1:30000][GF.df$'SH_GF_QC'[1:30000] == 3] <- NA
    GF.df$'SLE_FILLED'[1:30000][GF.df$'SLE_GF_QC'[1:30000] == 3] <- NA
    
    AEcombo_GF_df$'NR_G1_SH_SLE' <- AEcombo_GF_df$'NR_G1' - GF.df$'SH_FILLED' - GF.df$'SLE_FILLED'

    
    # . level_5: NETRAD + G1 + Sx + Spho
    AEcombo_GF_df$'NR_G1_SH_SLE_Spho' <- AEcombo_GF_df$'NR_G1_SH_SLE' - Spho
    
    # level_6: NETRAD + G1 + Sx + Spho + Sbio
    AEcombo_GF_df$'NR_G1_SH_SLE_Spho_Sbio' <- AEcombo_GF_df$'NR_G1_SH_SLE_Spho' - Sbio
    
    # Turbulent energy --------------------
    # Only 4 scenarios are possible given the availability of gapfilled data:
    # Hun_LEun, Hfk01_LEfk01, H_LE, Hust_LEust
    TEcombo_GF_df <- data.frame()
    
    # . level_0: UNCLEANED FLUXES
    
    # Rough filtering on ONEFlux GF FLAG == 3 
    # !! NOTA: si può fare perchè sono solo nella prima parte dell'anno, 
    GF.df$'H_UNCLEANED_FILLED'[1:30000][GF.df$'H_UNCLEANED_GF_QC'[1:30000] == 3] <- NA
    GF.df$'LE_UNCLEANED_FILLED'[1:30000][GF.df$'LE_UNCLEANED_GF_QC'[1:30000] == 3] <- NA
    
    TEcombo_GF_df <- data.frame(
      Hun_LEun = GF.df$'H_UNCLEANED_FILLED' + GF.df$'LE_UNCLEANED_FILLED'
    )

    if (any(TEcombo_GF_df$'Hun_LEun' > 5000, na.rm = TRUE)) {TEcombo_GF_df$'Hun_LEun'[TEcombo_GF_df$'Hun_LEun' > 5000] <- NA}
    if (any(TEcombo_GF_df$'Hun_LEun' < -5000, na.rm = TRUE)) {TEcombo_GF_df$'Hun_LEun'[TEcombo_GF_df$'Hun_LEun' < -5000] <- NA}
    if (any(TEcombo_GF_df$'Hun_LEun' > 5000, na.rm = TRUE) | any(TEcombo_GF_df$'Hun_LEun' < -5000, na.rm = TRUE)) {
      # re-fill it, but only internal gaps
      TEcombo_GF_df$'Hun_LEun' <- zoo::na.approx(TEcombo_GF_df$'Hun_LEun', na.rm=F) 
    }
    
    # . level_1: FLUXES (uncorrected) filtered by Foken2024 (FLAG == 2)
    if (ebc.ancillary$'network'[i.bd] == 'ICOS') {
      # Rough filtering on ONEFlux GF FLAG == 3
      GF.df$'H_fk01_FILLED'[1:30000][GF.df$'H_hf01_GF_QC'[1:30000] == 3] <- NA
      GF.df$'LE_UNCLEANED_FILLED'[1:30000][GF.df$'LE_UNCLEANED_GF_QC'[1:30000] == 3] <- NA
      
      TEcombo_GF_df$'Hfk01_LEfk01' <- GF.df$'H_fk01_FILLED' + GF.df$'LE_fk01_FILLED'
      
      if (any(TEcombo_GF_df$'Hfk01_LEfk01' > 3000, na.rm = TRUE)) {TEcombo_GF_df$'Hfk01_LEfk01'[TEcombo_GF_df$'Hfk01_LEfk01' > 3000] <- NA}
      if (any(TEcombo_GF_df$'Hfk01_LEfk01' < -3000, na.rm = TRUE)) {TEcombo_GF_df$'Hfk01_LEfk01'[TEcombo_GF_df$'Hfk01_LEfk01' < -3000] <- NA}
      if (any(TEcombo_GF_df$'Hfk01_LEfk01' > 3000, na.rm = TRUE) | any(TEcombo_GF_df$'Hfk01_LEfk01' < -3000, na.rm = TRUE)) {
        # re-fill it, but only internal gaps
        TEcombo_GF_df$'Hfk01_LEfk01' <- zoo::na.approx(TEcombo_GF_df$'Hfk01_LEfk01', na.rm=F) 
      }
    } else {
      
      TEcombo_GF_df$'Hfk01_LEfk01' <- NA
      
    }
    
    # . level_2: FULL CLEANED FLUXES
    
    # Rough filtering on ONEFlux GF FLAG == 3
    GF.df$'H_FILLED'[1:30000][GF.df$'H_GF_QC'[1:30000] == 3] <- NA
    GF.df$'LE_FILLED'[1:30000][GF.df$'LE_GF_QC'[1:30000] == 3] <- NA
    
    TEcombo_GF_df$'H_LE' <- GF.df$'H_FILLED' + GF.df$'LE_FILLED'
    
    if (any(TEcombo_GF_df$'H_LE' > 3000, na.rm = TRUE)) {TEcombo_GF_df$'H_LE'[TEcombo_GF_df$'H_LE' > 3000] <- NA}
    if (any(TEcombo_GF_df$'H_LE' < -3000, na.rm = TRUE)) {TEcombo_GF_df$'H_LE'[TEcombo_GF_df$'H_LE' < -3000] <- NA}
    if (any(TEcombo_GF_df$'H_LE' > 3000, na.rm = TRUE) | any(TEcombo_GF_df$'H_LE' < -3000, na.rm = TRUE)) {
      # re-fill it, but only internal gaps
      TEcombo_GF_df$'H_LE' <- zoo::na.approx(TEcombo_GF_df$'H_LE', na.rm=F) 
    }
    
    
    # . level_3: FULL CLEANED FLUXES USTAR FILTERED
    
    # Rough filtering on ONEFlux GF FLAG
    GF.df$'H_UF_FILLED'[1:30000][GF.df$'H_UF_GF_QC'[1:30000] == 3] <- NA
    GF.df$'LE_UF_FILLED'[1:30000][GF.df$'LE_UF_GF_QC'[1:30000] == 3] <- NA
    
    TEcombo_GF_df$'Hust_LEust' <- GF.df$'H_UF_FILLED' + GF.df$'LE_UF_FILLED'
    
    # Rough filtering on uncorrected fluxes sum
    if (any(TEcombo_GF_df$'Hust_LEust' > 5000, na.rm = TRUE)) {TEcombo_GF_df$'Hust_LEust'[TEcombo_GF_df$'Hust_LEust' > 5000] <- NA}
    if (any(TEcombo_GF_df$'Hust_LEust' < -5000, na.rm = TRUE)) {TEcombo_GF_df$'Hust_LEust'[TEcombo_GF_df$'Hust_LEust' < -5000] <- NA}
    if (any(TEcombo_GF_df$'Hust_LEust' > 5000, na.rm = TRUE) | any(TEcombo_GF_df$'Hust_LEust' < -5000, na.rm = TRUE)) {
      # re-fill it, but only internal gaps
      TEcombo_GF_df$'Hust_LEust' <- zoo::na.approx(TEcombo_GF_df$'Hust_LEust', na.rm=F) 
    }
    
  } else {
    
    AEcombo_GF_df <- data.frame()
    TEcombo_GF_df <- data.frame()
    cat(paste('\nNO gapfilled data for', cur.site, '\n'))
    
  }
  
  
  # EBC proxies ----------------------------------------------------------------
  # A dataframe with additional proxies to use as explaining variables for EBC dynamics
  
  ## Seasons ------------------------------------
  mesi <- as.numeric(format(strptime(base.data$'TIMESTAMP_START', '%Y%m%d%H%M', tz='GMT'), "%m"))
  seasons <- cut(mesi,
                 breaks = c(0, 2, 5, 8, 11, 12),
                 labels = c("DJF", "MAM", "JJA", "SON", "DJF"),
                 right = TRUE, include.lowest = TRUE)
  rm(mesi)
  
  ## GPP ------------------------------------
  if (!all(is.na(base.data$'GPP'))) {
    GPP.classes <- cut(base.data$'GPP',
                       breaks = c(0, seq(0.001, 50.001, 2.5)),
                       right = TRUE, include.lowest = TRUE)
    levels(GPP.classes) <- gsub(',','-', levels(GPP.classes)) # to avoid ',' in names 
  } else {
    GPP.classes <- rep(NA, nrow(base.data))
  }
  # table(GPP.classes)
  
  
  ## Stability  ---------------------------------
  # Atmospheric stratification by stability parameter
  # Atmospheric stability classes (Foken, 2008b)
  # Stability levels              Stability parameter
  # .................................................
  # very unstable                 -0.033 - -0.01
  # unstable                      -0.01 - -0.0033
  # neutral to slightly unstable  -0.0033 - 0.0002
  # neutral to slightly stable    0.0002 - 0.004
  # stable                        0.004 - 0.017
  # very stable                   > 0.017
  # vedi: Cui, W., Chui, T.F.M., 2019. Temporal and spatial variations of energy balance closure across FLUXNET research sites. Agric. For. Meteorol.
  ZL <- c(); ZL <- base.data$'ZL'
  atm.strat <- c()
  atm.strat <- cut(ZL,
                   breaks = c(min(ZL, na.rm = T), -0.01, -0.0033, 0.0002, 0.004, 0.017, max(ZL, na.rm = T)),
                   labels = c('very_unstable', 'unstable', 'neutral_unstable', 'neutral_stable', 'stable', 'very_stable'))
  rm(ZL)
  
  ## RH -----------------------------------------
  RH.classes <- c()
  RH.classes <- cut(replace(base.data$'RH', base.data$'RH' > 100, 100),
                    breaks = seq(0, 100, 5)
  )
  levels(RH.classes) <- gsub(',','-', levels(RH.classes)) # to avoid ',' in names

  
  ## Ustar --------------------------------------
  ust.tmp <- replace(base.data$'USTAR', base.data$'USTAR' < 0, 0)
  ust.tmp <- replace(base.data$'USTAR', base.data$'USTAR' > 3, 3)
  ust.classes <- c()
  ust.classes <- cut(ust.tmp,
                     breaks = seq(0, 3, 0.25))
  levels(ust.classes) <- gsub(',','-', levels(ust.classes)) # to avoid ',' in names
  rm(ust.tmp)
  
  ## Sun elevation & twilight -------------------
  # solar elevation hs = 90°− θz # θz = solar zenith
  # se hs > 0°→ Sole sopra l’orizzonte (giorno, geometricamente)
  # se hs ≤ 0° → Sole sotto l’orizzonte (notte, geometricamente)
  hs.classes <- cut(daynight_df$'elev_deg',
                    breaks = seq(-69, 69, 3) # Range completo teorico -90°, 90°, reale osservabile sulla Terra -70°,70°
  )
  levels(hs.classes) <- gsub(',','-', levels(hs.classes))  # to avoid ',' in names
  
  # Classi di crepuscolo:
  # Day: h > 0°
  # Alba/tramonto apparente: h = -0.833°
  # Civil twilight: −6° ≤ h ≤ 0°  (ancora luce diffusa)
  # Nautical twilight:  −12° ≤ h ≤ 6°   (orizzonte appena distinguibile)
  # Astronomical twilight:  −8° ≤ h ≤ 12°   (inizio/notte astronomica)
  # Night (astronomica):  h < 18°   (buio totale)
  twilight <- c()
  twilight <- daynight_df$'twilight'  
  
  
  # POXIES dataframe ..............
  base.proxies <- data.frame(NULL)
  base.proxies <- data.frame(
    
    # Seasons
    seasons = seasons,
    # GPP
    GPP_classes = GPP.classes,
    # Atmospheric stability conditions
    atm_strat = atm.strat,
    # RH
    RH_classes = RH.classes,
    # Ustar
    ustar_classes = ust.classes,
    # Sun Elevaton Angle
    SEA_classes = hs.classes,
    # Twilight
    twilight = twilight
  )
  
  # free some memory
  rm(seasons, GPP.classes, atm.strat, RH.classes, ust.classes, hs.classes, twilight)
  
  
  # EBC base data list ----------------------------------
  EBC.base.data[[i.bd]] <- vector(mode = 'list', length = 6)
  EBC.base.data[[i.bd]] <- list(
    TIMESTAMP_START = base.data$'TIMESTAMP_START',
    AE_combo = AEcombo_df,
    TE_combo = TEcombo_df,
    AE_GF_combo = AEcombo_GF_df,
    TE_GF_combo = TEcombo_GF_df,
    proxies = base.proxies
  )
  
  cat('\ndone!\n\n')
  
} # base data preparation loop END!

# save the "new data" 
new.data <- do.call(rbind, new.data)
fwrite(new.data, paste0(dirname(data.dir), "/newdata_Sbiopho_NRcorr_v", vernum, ".csv"), sep=',', quote = F, row.names = F)


# Save as Robjects
saveRDS(NETRAD.correction.summary, file = paste0(dirname(data.dir) , paste0("/NETRAD_corrections_stats_v", vernum, ".RData")))
saveRDS(EBC.base.data, file = paste0(dirname(data.dir) , paste0("/EBC_basedata_COMBOs_v", vernum, ".RData")))

gc()




# *** EBC COMPUTATION *** --------------------------------------

# Cross regressions, full-factorial, AE ~ TE
# output:
# RMA: slope, intercept, R², slope diff from 1
# OLS: slope, intercept, R², slope diff from 1
# OLS intercept = 0: slope, R², slope diff from 1
# EBC Imbalance 
# EBC Ratio (only gap-filled data) 


# eventualmente caricare i dati da qui:
if (!any(grepl('EBC.base.data', ls()))) { 
  EBC.base.data <- list()
  # EBC.base.data <- readRDS('EBC_DATA/EBC_basedata_COMBOs_v15.RData')  # previous
  EBC.base.data <- readRDS('EBC_DATA/EBC_basedata_COMBOs_v16.RData') # disponibile dopo aver verificato la tower disturbance
}

## EBC half-hour scale ----------------------
# . NON gap-filled data
# . perform:
#   . OLS
#   . OLS 0 intercept
#   . RMA 
#   . Imbalance 
#   . PROXIES analysis

EBC_results <- data.frame()

iebc <- 1
for (iebc in seq_along(EBC.base.data)) {   # loop over stations
  
  t0 <- Sys.time()
  
  # Estrai i dati per la stazione corrente
  cur.site <- names(EBC.base.data)[iebc]
  cur.ntwrk <- ebc.ancillary$'network'[ebc.ancillary$'station' == cur.site]
  cur.igbp <- ebc.ancillary$'IGBP'[ebc.ancillary$'station' == cur.site]
  
  cat(paste0('\n+++++++++++++ Computing EBC for ', cur.site, ' +++++++++++++\n'))
  
  AE <- EBC.base.data[[cur.site]]$'AE_combo'
  TE <- EBC.base.data[[cur.site]]$'TE_combo'
  
  
  ### ---- Scala semioraria full dataset ----
  
  cat('Scala semioraria full dataset:\n')
  
  for (ae_var in names(AE)) {
    
    for (te_var in names(TE)) {
      
      dati <- data.frame()
      dati <- data.frame(x = AE[[ae_var]], y = TE[[te_var]])
      dati <- na.omit(dati)
      
      if (nrow(dati) > 2) {
        
        # EBC computation ......................................................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance: IMB=NR-AE-TE
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        ebc.imb <- mean(AE[[ae_var]] - TE[[te_var]], na.rm=T) # mean
        ebc.imb.se <- sd(AE[[ae_var]] - TE[[te_var]], na.rm=T) / sqrt(length(na.omit(AE[[ae_var]] - TE[[te_var]]))) # SE
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE)*100,2)
        
        # Store results
        EBC_results <- bind_rows(EBC_results, 
                                 cbind(
                                   site = cur.site,
                                   network = cur.ntwrk,
                                   igbp = cur.igbp,
                                   t_scale = "semioraria",
                                   fattore = NA,
                                   livello = NA,
                                   AE_var = ae_var,
                                   TE_var = te_var,
                                   res,
                                   EBC_imb = ebc.imb,
                                   EBC_imb_se = ebc.imb.se,
                                   datacov_a = datacov_a,
                                   datacov_p = datacov_p
                                 )
        )
      } else {
        # print(paste0('NO DATA FOR: ', ae_var,'-',te_var))
      }
      
    } # regression loop over TE END
    
  } # regression loop over AE END
  
  cat('done.\n')
  
  
  ## ---- Per proxies ----
  # full dataset
  
  cat('Scala semioraria per proxies:\n')
  
  # Subset the proxies dataset (retain only factors, which want to analyse)
  # SELECT OTHERS HERE IF NECESSARY
  proxies <- data.frame()
  proxies <- EBC.base.data[[cur.site]]$'proxies'
  proxies <- proxies[, names(proxies) %in% c("GPP_classes","RH_classes","atm_strat","SEA_classes","twilight")]       # names(proxies)
  
  # Select the AE and TE combos to use for the proxies analysis
  # ADD OTHERS HERE IF NECESSARY
  AE.proxy.sel <- c("NR","NR_G1_SH_SLE_Spho_Sbio", "NRslope_G1_SH_SLE_Spho_Sbio")    # names(EBC.base.data[[iebc]]$AE_combo)
  TE.proxy.sel <- c("H_LE","Hspanal_LEspanal","Hust_LEust")                                        # names(EBC.base.data[[iebc]]$TE_combo)
  
  for (fattore in names(proxies)) {   # regression loop over proxies
    
    for (liv in unique(proxies[[fattore]])) {  # regression loop over proxie's levels
      
      idx <- which(proxies[[fattore]] == liv)
      
      if (length(idx) > 2) {
        
        for (ae_var in AE.proxy.sel) {  # regression loop over AE
          
          for (te_var in TE.proxy.sel) {  # regression loop over TE 
            
            dati <- data.frame(x = AE[[ae_var]][idx], y = TE[[te_var]][idx]) # this is possible as AE, TE and proxies have the same length
            dati <- na.omit(dati)
            
            if (nrow(dati) > 2) {
              
              # EBC computation ................................................
              
              # 1. Regressions: RMA and OLS (free and 0 intercept))
              res <- data.frame()
              res <- calc_reg(dati$x, dati$y) 
              
              # 2. EBC Imbalance 
              ebc.imb <- NULL; ebc.imb.se <- NULL 
              ebc.imb <- mean(AE[[ae_var]][idx] - TE[[te_var]][idx], na.rm=T) # mean
              ebc.imb.se <- sd(AE[[ae_var]][idx] - TE[[te_var]][idx], na.rm=T) / sqrt(length(na.omit(AE[[ae_var]][idx] - TE[[te_var]][idx]))) # SE
              
              # Data coverage
              datacov_a <- nrow(dati)
              datacov_p <- round(nrow(dati)/nrow(AE)*100,2)
              
              # Store results
              EBC_results <- bind_rows(EBC_results, 
                                       cbind(
                                         site = cur.site,
                                         network = cur.ntwrk,
                                         igbp = cur.igbp,
                                         t_scale = "semioraria",
                                         fattore = fattore,
                                         livello = liv,
                                         AE_var = ae_var,
                                         TE_var = te_var,
                                         res,
                                         EBC_imb = ebc.imb,
                                         EBC_imb_se = ebc.imb.se,
                                         datacov_a = datacov_a,
                                         datacov_p = datacov_p
                                       ))
            }
          }  # regression loop over TE END
          
        }  # regression loop over AE END
        
      }
    }   # regression loop over proxie's levels END
    
  }   # regression loop over proxies END
  
  cat('done.\n')
  
  
  ## ---- 2024 only ----
  # needed mostly to check the consistency with the full-period datasets, and to analyse the NR corrections 
  # (for consistency with the FoV correction as it is implemented for 2024 only)
  
  cat('2024:\n')
  
  # Trim the AE eand TE dataset to the desired period
  t1 <- which(EBC.base.data[[cur.site]]$'TIMESTAMP_START' == '202401010000')
  t2 <- which(EBC.base.data[[cur.site]]$'TIMESTAMP_START' == '202412312330')
  if (length(t2) == 0) {t2 = nrow(AE)}
  AE_2024 <- data.frame()
  AE_2024 <- AE[t1:t2, ]
  TE_2024 <- data.frame()
  TE_2024 <- TE[t1:t2, ]
  
  for (ae_var in names(AE_2024)) {
    
    for (te_var in names(TE_2024)) {
      
      dati <- data.frame()
      dati <- data.frame(x = AE_2024[[ae_var]], y = TE_2024[[te_var]])
      dati <- na.omit(dati)
      
      if (nrow(dati) > 2) {
        
        # EBC computation ......................................................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        ebc.imb <- mean(AE_2024[[ae_var]] - TE_2024[[te_var]], na.rm=T) # mean
        ebc.imb.se <- sd(AE_2024[[ae_var]] - TE_2024[[te_var]], na.rm=T) / sqrt(length(na.omit(AE_2024[[ae_var]] - TE_2024[[te_var]]))) # SE
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE_2024)*100,2)
        
        EBC_results <- bind_rows(EBC_results, 
                                 cbind(
                                   site = cur.site,
                                   network = cur.ntwrk,
                                   igbp = cur.igbp,
                                   t_scale = "2024",
                                   fattore = NA,
                                   livello = NA,
                                   AE_var = ae_var,
                                   TE_var = te_var,
                                   res,
                                   EBC_imb = ebc.imb,
                                   EBC_imb_se = ebc.imb.se,
                                   datacov_a = datacov_a,
                                   datacov_p = datacov_p
                                 )
        )
      } else {
        # print(paste0('\nNO DATA FOR: ', ae_var,'-',te_var))
      }
    }
  } # 2024 loop end
  
  rm(AE_2024, TE_2024)
  
  
  cat('done.\n')
  
  tx <- Sys.time()
  ellapsed.time <- difftime(tx, t0, units="hours")
  cat(paste('\n......', cur.site, 'EBC regressions computed in:', round(ellapsed.time, 3), 'hours ......\n\n'))
  
}    # loop over stations END

# Save it!
fwrite(EBC_results, paste0('OUTPUT/EBC_results_HH_nofill_v', vernum, '.csv'), sep=',', na = NA, quote = F, row.names = F)




## EBC multiple time scale -------------------------------
# . Gap-filled data
# . perform OLS, OLS 0 intercept, RMA, EBCratio, (EBCabs?), and temporal aggregations analysis

# Associa a ogni data la stagione 
season_map <- function(dates) {
  m <- lubridate::month(dates)
  s <- case_when(
    m %in% c(12, 1, 2) ~ "DJF",
    m %in% c(3, 4, 5) ~ "MAM",
    m %in% c(6, 7, 8) ~ "JJA",
    m %in% c(9, 10, 11) ~ "SON"
  )
  return(factor(s, levels = c("DJF", "MAM", "JJA", "SON")))
}

# Set how to compute the energy imbalance: TRUE == mean, FALSE == median
imb.as.mean <- FALSE

# hosting data structure
EBC_results_GF <- data.frame()

iebc <- 1
for (iebc in seq_along(EBC.base.data)) {   # loop over stations
  
  t0 <- Sys.time()
  
  # Estrai i dati per la stazione corrente
  cur.site <- names(EBC.base.data)[iebc]
  cur.ntwrk <- ebc.ancillary$'network'[ebc.ancillary$'station' == cur.site]
  cur.igbp <- ebc.ancillary$'IGBP'[ebc.ancillary$'station' == cur.site]
  
  cat(paste0('\n+++++++++++++ Computing EBC (gap-filled data) for ', cur.site, ' +++++++++++++\n'))
  
  AE <- data.frame()
  AE <- EBC.base.data[[cur.site]]$'AE_GF_combo'
  TE <- data.frame()
  TE <- EBC.base.data[[cur.site]]$'TE_GF_combo'
  
  # NOTE_202508: Proxies analysis is NOT performed over GF data (for the moment)
  
  
  ### ---- Scala semioraria  ----
  # EBC stats calcolate a scala semioraria, su tutto il dataset
  # - regressioni 
  # - EBC Imbalance
  # - EBC ratio (bulk EBC fraction)
  
  cat('Scala semioraria (gap-filled):\n')
  
  for (ae_var in names(AE)) {  # regression loop over AE
    
    for (te_var in names(TE)) {  # regression loop over TE
      
      dati <- data.frame()
      dati <- data.frame(x = AE[[ae_var]], y = TE[[te_var]])
      dati <- na.omit(dati)
      
      if (nrow(dati) > 2 & !any(apply(dati, 2, sum, na.rm=T) == 0)) { # This is mainly for NEON data that does not have foken filtered data
        
        # EBC computation ...........................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        if (imb.as.mean) {
          ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
        } else {
          ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
        }  
        ebc.imb.se <- sd(AE[[ae_var]] - TE[[te_var]], na.rm=T) / sqrt(length(na.omit(AE[[ae_var]] - TE[[te_var]]))) # SE
        
        # 3. EBC ratio: TE/AE
        ebc.ratio <- NULL
        ebc.ratio <- sum(dati$y, na.rm=TRUE) / sum(dati$x, na.rm=T)
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE)*100,2)
        
        # Fill the output
        EBC_results_GF <- bind_rows(EBC_results_GF, 
                                    cbind(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "semioraria",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var,
                                      res,
                                      EBC_imb = ebc.imb,
                                      EBC_imb_se = ebc.imb.se,
                                      EBC_imb_abs = NA,
                                      EBC_imb_abs_se = NA,
                                      EBC_imb_cum = NA,
                                      EBC_imb_cum_ciL = NA,
                                      EBC_imb_cum_ciU = NA,
                                      EBC_ratio = ebc.ratio,
                                      datacov_a = datacov_a,
                                      datacov_p = datacov_p
                                    )
        )
        
      } else {
        
        EBC_results_GF <- bind_rows(EBC_results_GF, 
                                    data.frame(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "semioraria",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var
                                    )
        )
        
      }
    }
    
  } # Half-hourly data loop end
  
  cat('done.\n')
  
  
  # ---- Scala semioraria/stagione ----
  # EBC stats calcolate a scala semioraria, per stagione
  # - regressioni 
  # - EBC Imbalance
  # - EBC ratio (bulk EBC fraction)
  
  cat('Scala semioraria per stagione (gap-filled):\n')
  
  # Timestamp in formato POSIXct (TS.cur)
  TS.cur <- lubridate::ymd_hm(as.character(EBC.base.data[[iebc]]$'TIMESTAMP_START'))
  
  # Crea vettore stagioni per ciascun semiorario
  TS.season <- season_map(TS.cur)
  
  # AE stagionali (assegnazione della stagione ad ogni semiora)
  AE_season <- AE %>%
    mutate(season = TS.season) %>%
    # group_by(season) %>%
    # summarise(across(everything(), sum, na.rm = TRUE)) %>%
    # ungroup() %>%
    as.data.frame()
  
  # TE stagionali (assegnazione della stagione ad ogni semiora)
  TE_season <- TE %>%
    mutate(season = TS.season) %>%
    # group_by(season) %>%
    # summarise(across(everything(), sum, na.rm = TRUE)) %>%
    # ungroup() %>%
    as.data.frame()
  
  
  # Loop per regressioni stagionali
  for (season in levels(TS.season)) {
    
    for (ae_var in names(AE)) {
      
      for (te_var in names(TE)) {
        
        # EBC computation ...........................
                
        # Dati stagionali per le due variabili
        dati <- data.frame(
          x = AE_season[AE_season$season == season, ae_var],
          y = TE_season[TE_season$season == season, te_var]
        )
        dati <- na.omit(dati)
        
        if (nrow(dati) > 2 & !any(apply(dati, 2, sum, na.rm = TRUE) == 0)) {
          
          # 1. Regressioni RMA + OLS
          res <- data.frame()
          res <- calc_reg(dati$x, dati$y)
          
          # 2. EBC Imbalance
          ebc.imb <- NULL; ebc.imb.se <- NULL 
          if (imb.as.mean) {
            ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
          } else {
            ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
          }
          ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
          
          # 3. EBC ratio (stagionale)
          ebc.ratio <- sum(dati$y, na.rm = TRUE) / sum(dati$x, na.rm = TRUE)
          
          # 4. Data coverage
          datacov_a <- nrow(dati)
          datacov_p <- round(nrow(dati) / nrow(AE_season[AE_season$season == season, ]) * 100, 2)
          
          # Aggiungi riga a EBC_results_GF
          EBC_results_GF <- bind_rows(EBC_results_GF,
                                      cbind(
                                        site = cur.site,
                                        network = cur.ntwrk,
                                        igbp = cur.igbp,
                                        t_scale = "semioraria_stagionale",
                                        fattore = "season",
                                        livello = season,
                                        AE_var = ae_var,
                                        TE_var = te_var,
                                        res,
                                        EBC_imb = ebc.imb,
                                        EBC_imb_se = ebc.imb.se,
                                        EBC_imb_abs = NA,
                                        EBC_imb_abs_se = NA,
                                        EBC_imb_cum = NA,
                                        EBC_imb_cum_ciL = NA,
                                        EBC_imb_cum_ciU = NA,
                                        EBC_ratio = ebc.ratio,
                                        datacov_a = datacov_a,
                                        datacov_p = datacov_p
                                      ))
          
        } else {
          EBC_results_GF <- bind_rows(EBC_results_GF,
                                      data.frame(
                                        site = cur.site,
                                        network = cur.ntwrk,
                                        igbp = cur.igbp,
                                        t_scale = "semioraria_stagionale",
                                        fattore = "season",
                                        livello = NA,
                                        AE_var = ae_var,
                                        TE_var = te_var
                                      ))
        }
      }
      
    } # loop semiore end  
    
  } # loop stagione end
  
  
  # ---- Scala semioraria/stagione/day-night ----
  # EBC stats calcolate a scala semioraria, per stagione e per daylight
  # NOTE: this will be used for the IMB vs GPP analysis
  
  cat('Scala semioraria per stagione per daylight (gap-filled):\n')
    
  # Converto il proxi twilight in fattore
  twilight <- as.factor(EBC.base.data[[cur.site]]$'proxies'$'twilight')
  
  # check
  stopifnot(length(twilight) == length(TS.cur))
  
  # AE e TE con stagione + twilight
  AE_season_dl <- data.frame()
  AE_season_dl <- AE %>%
    mutate(
      season   = TS.season,
      twilight = twilight,
      year = lubridate::year(TS.cur) # colonna con l'anno di riferimento, serve per IMB cumulato
    ) %>%
    as.data.frame()
  
  TE_season_dl <- TE %>%
    mutate(
      season   = TS.season,
      twilight = twilight
    ) %>%
    as.data.frame()
  
  
  # Loop per regressioni per combinazione stagione x twilight
  for (season in levels(TS.season)) {
    # season = levels(TS.season)[1]
    
    for (tw in levels(twilight)) {
      # tw = levels(twilight)[1]
      
      for (ae_var in names(AE)) {
        # ae_var = names(AE)[7]
        
        for (te_var in names(TE)) {
          # te_var = names(TE)[3]
          
          # Filtra i dati per stagione + twilight
          dati <- data.frame(
            x = AE_season_dl[AE_season_dl$season == season & AE_season_dl$twilight == tw, ae_var],
            y = TE_season_dl[TE_season_dl$season == season & TE_season_dl$twilight == tw, te_var],
            year = AE_season_dl[AE_season_dl$season == season & AE_season_dl$twilight == tw, 'year']
          )
          colnames(dati) <- c("x","y","year")
          dati <- na.omit(dati)
          
          if (nrow(dati) > 2 && !any(apply(dati, 2, sum, na.rm = TRUE) == 0)) {
            
            # 1. Regressioni RMA + OLS
            res <- calc_reg(dati$x, dati$y)
            
            # 2. EBC Imbalance (mean o median)
            if (imb.as.mean) {
              ebc.imb <- mean(dati$x - dati$y, na.rm = TRUE)
            } else {
              ebc.imb <- median(dati$x - dati$y, na.rm = TRUE)
            }
            ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) /
              sqrt(length(na.omit(dati$x - dati$y)))
            
            # 2.1. CUMULATIVE EBC Imbalance (mean)
            # NOTE_20251128: solo per questa analisi (IMB ~ GPP) l'imbalance è calcolato come CUMULATA
            #                OUTPUT: IMB medio delle cumulate stagionali + twilight di ogni anno ± CI 95% bootstrapped
            # Cumulative IMB per year/season/twilight
            ebc.imb_cum.y <- tapply(dati$x - dati$y, dati$year, sum)
            # # Mean cumulative IMB per year/season/twilight: DATI REALI

            # Bootstrapping
            set.seed(123)
            B <- 2000 # bootstrap replicates
            # bootstrap per la media
            boot_means  <- replicate(B, mean(sample(ebc.imb_cum.y, replace = TRUE)))
            # Mean cumulative IMB per year/season/twilight: DATI BOOTs
            ebc.imb_cum <- mean(boot_means)
            # CI 95%
            ci_lower <- quantile(boot_means, c(0.025))
            ci_upper <- quantile(boot_means, c(0.975))
            
            # 3. EBC ratio
            ebc.ratio <- sum(dati$y, na.rm = TRUE) / sum(dati$x, na.rm = TRUE)
            
            # 4. Data coverage (rispetto a tutte le semiore della season + twilight)
            n_tot <- nrow(AE_season_dl[AE_season_dl$season == season &
                                         AE_season_dl$twilight == tw, ])
            datacov_a <- nrow(dati)
            datacov_p <- if (n_tot > 0) round(datacov_a / n_tot * 100, 2) else NA
            
            # Aggiungi riga a EBC_results_GF
            EBC_results_GF <- bind_rows(
              EBC_results_GF,
              cbind(
                site   = cur.site,
                network = cur.ntwrk,
                igbp   = cur.igbp,
                t_scale = "semioraria_stagionale_dl",
                fattore = "season_dl",
                livello = paste0(season, '_', tw),     # livello = season_twilight
                AE_var = ae_var,
                TE_var = te_var,
                res,
                EBC_imb      = ebc.imb,
                EBC_imb_se   = ebc.imb.se,
                EBC_imb_abs  = NA,
                EBC_imb_abs_se = NA,
                EBC_imb_cum = ebc.imb_cum,
                EBC_imb_cum_ciL = ci_lower,
                EBC_imb_cum_ciU = ci_upper,
                EBC_ratio    = ebc.ratio,
                datacov_a    = datacov_a,
                datacov_p    = datacov_p
              )
            )
            
          } else {
            
            EBC_results_GF <- bind_rows(
              EBC_results_GF,
              data.frame(
                site    = cur.site,
                network = cur.ntwrk,
                igbp    = cur.igbp,
                t_scale = "semioraria_stagionale_dl",
                fattore = "season_dl",
                livello = paste0(season, '_', tw),
                AE_var  = ae_var,
                TE_var  = te_var
              )
            )
          }
          
        } # loop te_var
        
      } # loop ae_var
      
    } # loop twilight
    
  } # loop season
  
  
  
  ### ---- Scala semioraria/annuale ----
  # EBC stats calcolate a scala semioraria, per anno
  # NOTE: this will be used for the IMB vs BIOMASS analysis
  
  cat('Scala semioraria per anno (gap-filled):\n')
  
  # # Timestamp in formato POSIXct (TS.cur)
  # TS.cur <- lubridate::ymd_hm(as.character(EBC.base.data[[iebc]]$'TIMESTAMP_START'))
  
  # Crea vettore anno per ciascun semiorario
  TS.year <- c()
  TS.year <- lubridate::year(TS.cur)
  
  # AE per anno (assegnazione dell' anno ad ogni semiora)
  AE_year <- AE %>%
    mutate(year = TS.year) %>%
    # group_by(year) %>%
    # summarise(across(everything(), sum, na.rm = TRUE)) %>%
    # ungroup() %>%
    as.data.frame()
  
  # TE per anno (assegnazione dell' anno ad ogni semiora)
  TE_year <- TE %>%
    mutate(year = TS.year) %>%
    # group_by(year) %>%
    # summarise(across(everything(), sum, na.rm = TRUE)) %>%
    # ungroup() %>%
    as.data.frame()
  
  
  # Loop per regressioni annuali
  for (year in unique(TS.year)) {
    
    for (ae_var in names(AE)) {
      
      for (te_var in names(TE)) {
        
        # EBC computation ...........................
        
        # Dati stagionali per le due variabili
        dati <- data.frame(
          x = AE_year[AE_year$year == year, ae_var],
          y = TE_year[TE_year$year == year, te_var]
        )
        dati <- na.omit(dati)
        
        if (nrow(dati) > 2 & !any(apply(dati, 2, sum, na.rm = TRUE) == 0)) {
          
          # 1. Regressioni RMA + OLS
          res <- data.frame()
          res <- calc_reg(dati$x, dati$y)
          
          # 2. EBC Imbalance
          ebc.imb <- NULL; ebc.imb.se <- NULL 
          if (imb.as.mean) {
            ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
          } else {
            ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
          }
          ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
          
          # 2.1. ABSOLUTE EBC Imbalance 
          # NOTE: this will be used for the IMB ~ BIOMASS analysis
          
          ebc.imb_abs <- NULL; ebc.imb.se_abs <- NULL 
          if (imb.as.mean) {
            ebc.imb_abs <- mean(abs(dati$x - dati$y), na.rm=T) # mean
          } else {
            ebc.imb_abs <- median(abs(dati$x - dati$y), na.rm=T) # median
          }
          ebc.imb.se_abs <- sd(abs(dati$x - dati$y), na.rm = TRUE) / sqrt(length(na.omit(abs(dati$x - dati$y))))
          
          # 3. EBC ratio (stagionale)
          ebc.ratio <- sum(dati$y, na.rm = TRUE) / sum(dati$x, na.rm = TRUE)
          
          # 4. Data coverage
          datacov_a <- nrow(dati)
          datacov_p <- round(nrow(dati) / nrow(AE_year[AE_year$year == year, ]) * 100, 2)
          
          # Aggiungi riga a EBC_results_GF
          EBC_results_GF <- bind_rows(EBC_results_GF,
                                      cbind(
                                        site = cur.site,
                                        network = cur.ntwrk,
                                        igbp = cur.igbp,
                                        t_scale = "semioraria_annuale",
                                        fattore = "year",
                                        livello = as.character(year),
                                        AE_var = ae_var,
                                        TE_var = te_var,
                                        res,
                                        EBC_imb = ebc.imb,
                                        EBC_imb_se = ebc.imb.se,
                                        EBC_imb_abs = ebc.imb_abs,
                                        EBC_imb_abs_se = ebc.imb.se_abs,
                                        EBC_imb_cum = NA,
                                        EBC_imb_cum_ciL = NA,
                                        EBC_imb_cum_ciU = NA,
                                        EBC_ratio = ebc.ratio,
                                        datacov_a = datacov_a,
                                        datacov_p = datacov_p
                                      ))
          
        } else {
          EBC_results_GF <- bind_rows(EBC_results_GF,
                                      data.frame(
                                        site = cur.site,
                                        network = cur.ntwrk,
                                        igbp = cur.igbp,
                                        t_scale = "semioraria_annuale",
                                        fattore = "year",
                                        livello = NA,
                                        AE_var = ae_var,
                                        TE_var = te_var
                                      ))
        }
      }
      
    } # loop semiore end  
    
  } # loop stagione end
  
  
  ## ---- Scala giornaliera ----
  # EBC stats calcolate a scala giornaliera, VALORI AGGRAGATI
  
  cat('Scala giornaliera (gap-filled):\n')
  
  # Converte timestamp in POSIXct
  # TS.cur <- c()
  # TS.cur <- EBC.base.data[[iebc]]$'TIMESTAMP_START'
  # TS.cur <- lubridate::ymd_hm(as.character(TS.cur))
  
  # Aggregazione giornaliere di AE ..........
  TS.day <- lubridate::as_date(TS.cur) # formattazione come giorno fatta qui con as.Date
  AE_day <- data.frame()
  AE_day <- AE %>% 
    group_by(day = TS.day) %>%                          
    summarise(across(everything(), sum, na.rm = TRUE))  # DAILY SUMS
  AE_day <- as.data.frame(AE_day)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  AE_day <- AE_day %>% 
    mutate_if(is.numeric, ~na_if(., 0))

  
  
  # Aggregazione giornaliere di TE ..........
  TE_day <- data.frame()
  TE_day <- TE %>%
    group_by(day = TS.day) %>%
    summarise(across(everything(), sum, na.rm = TRUE))
  TE_day <- as.data.frame(TE_day)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  TE_day <- TE_day %>% mutate_if(is.numeric, ~na_if(., 0))
  
  
  for (ae_var in names(AE)) { # regression loop over AE
    
    for (te_var in names(TE)) { # regression loop over TE 
      
      # ae_var <- names(AE)[1]; te_var <- names(TE)[1];
      dati <- data.frame(
        x = AE_day[[ae_var]],
        y = TE_day[[te_var]]
      )
      dati <- na.omit(dati)
      
      if (nrow(dati) > 2 & !any(apply(dati, 2, sum, na.rm=T) == 0)) { # This is mainly for NEON data that does not have foken filtered data
        
        # EBC computation ...........................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        if (imb.as.mean) {
          ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
        } else {
          ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
        }  
        ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
        
        # 3. EBC ratio: TE/AE
        ebc.ratio <- NULL
        ebc.ratio <- sum(dati$y, na.rm=TRUE) / sum(dati$x, na.rm=T)
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE_day)*100,2)
        
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    cbind(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "giornaliera",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var,
                                      res,
                                      EBC_imb = ebc.imb,
                                      EBC_imb_se = ebc.imb.se,
                                      EBC_imb_abs = NA,
                                      EBC_imb_abs_se = NA,
                                      EBC_imb_cum = NA,
                                      EBC_imb_cum_ciL = NA,
                                      EBC_imb_cum_ciU = NA,
                                      EBC_ratio = ebc.ratio,
                                      datacov_a = datacov_a,
                                      datacov_p = datacov_p
                                    )
        )
        
      } else { 
        
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    data.frame(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "giornaliera",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var
                                    )
        )
      }
      
    } # regression loop over TE END
    
  } # regression loop over AE END
  
  cat('done.\n')
  
  
  ## ---- Scala settimanale ----
  # EBC stats calcolate a scala settimanale, VALORI AGGREGATI
  cat('Scala settimanale (gap-filled):\n')
  
  # Converte timestamp in POSIXct
  # TS.cur <- c()
  # TS.cur <- EBC.base.data[[iebc]]$'TIMESTAMP_START'
  # TS.cur <- lubridate::ymd_hm(as.character(TS.cur))
  
  # Aggregazione settimanale di AE
  # Format timestamps as week of the year
  TS.week <- format(TS.cur, '%Y%W') 
  AE_week <- data.frame()
  AE_week <- AE %>%
    group_by(week = TS.week) %>%
    summarise(across(everything(), sum, na.rm = TRUE)) # WEEKLY SUMS
  AE_week <- as.data.frame(AE_week)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  AE_week <- AE_week %>% mutate_if(is.numeric, ~na_if(., 0))
  
  # Aggregazione settimanale di TE
  TE_week <- data.frame()
  TE_week <- TE %>%
    group_by(week = TS.week) %>%
    summarise(across(everything(), sum, na.rm = TRUE))
  TE_week <- as.data.frame(TE_week)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  TE_week <- TE_week %>% mutate_if(is.numeric, ~na_if(., 0))
  
  
  for (ae_var in names(AE)) { # regression loop over AE
    
    for (te_var in names(TE)) { # regression loop over TE 
      
      # ae_var <- names(AE)[1]; te_var <- names(TE)[1];
      dati <- data.frame(NULL)
      dati <- data.frame(
        x = AE_week[[ae_var]],
        y = TE_week[[te_var]]
      )
      dati <- na.omit(dati)
      
      if (nrow(dati) > 2 & !any(apply(dati, 2, sum, na.rm=T) == 0)) {
        
        # EBC computation ......................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        if (imb.as.mean) {
          ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
        } else {
          ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
        }  
        ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
        
        # 3. EBC ratio: TE/AE
        ebc.ratio <- NULL
        ebc.ratio <- sum(dati$y, na.rm=TRUE) / sum(dati$x, na.rm=T)
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE_week)*100,2)
        
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    cbind(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "settimanale",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var,
                                      res,
                                      EBC_imb = ebc.imb,
                                      EBC_imb_se = ebc.imb.se,
                                      EBC_imb_abs = NA,
                                      EBC_imb_abs_se = NA,
                                      EBC_imb_cum = NA,
                                      EBC_imb_cum_ciL = NA,
                                      EBC_imb_cum_ciU = NA,
                                      EBC_ratio = ebc.ratio,
                                      datacov_a = datacov_a,
                                      datacov_p = datacov_p
                                    )
        )
      } else {
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    data.frame(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "settimanale",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var
                                    )
        )
      }
      
    } # regression loop over TE END
    
  } # regression loop over AE END
    
  cat('done.\n')
  
  
  # ---- Scala mensile  ----
  # EBC stats calcolate a scala mensile, VALORI AGGRAGATI
  
  cat('Scala mensile (gap-filled):\n')
  
  # Converte timestamp in POSIXct
  # TS.cur <- c()
  # TS.cur <- EBC.base.data[[iebc]]$'TIMESTAMP_START'
  # TS.cur <- lubridate::ymd_hm(as.character(TS.cur))
  
  # Aggregazione mensile di AE
  # Format timestamps as months
  TS.mon <- format(TS.cur, '%Y%m') # CHECKIT OUT!
  AE_month <- data.frame()
  AE_month <- AE %>%
    group_by(month = TS.mon) %>%
    summarise(across(everything(), sum, na.rm = TRUE))  # MONTHLY SUMS
  AE_month <- as.data.frame(AE_month)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  AE_month <- AE_month %>% mutate_if(is.numeric, ~na_if(., 0))
  
  # Aggregazione settimanale di TE
  TE_month <- data.frame()
  TE_month <- TE %>%
    group_by(month = TS.mon) %>%
    summarise(across(everything(), sum, na.rm = TRUE))
  TE_month <- as.data.frame(TE_month)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  TE_month <- TE_month %>% mutate_if(is.numeric, ~na_if(., 0))
  
  
  for (ae_var in names(AE)) { # regression loop over AE
    
    for (te_var in names(TE)) { # regression loop over TE 
      
      dati <- data.frame(NULL)
      dati <- data.frame(
        x = AE_month[[ae_var]],
        y = TE_month[[te_var]]
      )
      dati <- na.omit(dati)
      
      if (nrow(dati) > 2 & !any(apply(dati, 2, sum, na.rm=T) == 0)) {
        
        # EBC computation ......................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        if (imb.as.mean) {
          ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
        } else {
          ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
        }  
        ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
        
        # 3. EBC ratio: TE/AE
        ebc.ratio <- NULL
        ebc.ratio <- sum(dati$y, na.rm=TRUE) / sum(dati$x, na.rm=T)
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE_month)*100,2)
        
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    cbind(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "mensile",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var,
                                      res,
                                      EBC_imb = ebc.imb,
                                      EBC_imb_se = ebc.imb.se,
                                      EBC_imb_abs = NA,
                                      EBC_imb_abs_se = NA,
                                      EBC_imb_cum = NA,
                                      EBC_imb_cum_ciL = NA,
                                      EBC_imb_cum_ciU = NA,
                                      EBC_ratio = ebc.ratio,
                                      datacov_a = datacov_a,
                                      datacov_p = datacov_p
                                    )
        )
      } else {
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    data.frame(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "mensile",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var
                                    )
        )
      }
      
    } # regression loop over TE END
    
  } # regression loop over AE END
  
  cat('done.\n')
  
  
  ## ---- Scala stagionale  ----
  # EBC stats calcolate a scala stagionale, VALORI AGGRAGATI
  
  cat('Scala stagionale (gap-filled):\n')
  
  # # Converte timestamp in POSIXct
  # TS.cur <- c()
  # TS.cur <- EBC.base.data[[iebc]]$'TIMESTAMP_START'
  # TS.cur <- lubridate::ymd_hm(as.character(TS.cur))
  
  # Aggregazione stagionale di AE
  AE_season <- data.frame()
  AE_season <- AE %>%
    mutate(
      month = month(TS.cur),
      # define HERE the season's months 
      season = case_when(
        month %>% str_detect("12|1|2") ~ "winter",
        month %>% str_detect("3|4|5") ~ "spring",
        month %>% str_detect("6|7|8") ~ "summer",
        TRUE ~ "autumn"
      )
    ) %>%
    group_by(season = season) %>%
    summarise(across(everything(), sum, na.rm = TRUE))  # SEASONAL SUMS
  AE_season <- as.data.frame(AE_season)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  AE_season <- AE_season %>% mutate_if(is.numeric, ~na_if(., 0))
  
  # Aggregazione stagionale di AE
  TE_season <- data.frame()
  TE_season <- TE %>%
    mutate(
      month = month(TS.cur),
      # define HERE the season's months 
      season = case_when(
        month %>% str_detect("12|1|2") ~ "winter",
        month %>% str_detect("3|4|5") ~ "spring",
        month %>% str_detect("6|7|8") ~ "summer",
        TRUE ~ "autumn"
      )
    ) %>%
    group_by(season = season) %>%
    summarise(across(everything(), sum, na.rm = TRUE))
  TE_season <- as.data.frame(TE_season)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  TE_season <- TE_season %>% mutate_if(is.numeric, ~na_if(., 0))
  
  
  for (ae_var in names(AE)) { # regression loop over AE
    
    for (te_var in names(TE)) { # regression loop over TE 
      
      # ae_var <- names(AE)[1]; te_var <- names(TE)[1];
      dati <- data.frame(NULL)
      dati <- data.frame(
        x = AE_season[[ae_var]],
        y = TE_season[[te_var]]
      )
      dati <- na.omit(dati)
      
      if (nrow(dati) > 0 & !any(apply(dati, 2, sum, na.rm=T) == 0)) {
        
        # EBC computation ......................................................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        res <- data.frame()
        res <- calc_reg(dati$x, dati$y) 
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        if (imb.as.mean) {
          ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
        } else {
          ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
        }  
        ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
        
        # 3. EBC ratio: TE/AE
        ebc.ratio <- NULL
        ebc.ratio <- sum(dati$y, na.rm=TRUE) / sum(dati$x, na.rm=T)
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE_season)*100,2)
        
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    cbind(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "stagionale",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var,
                                      res,
                                      EBC_imb = ebc.imb,
                                      EBC_imb_se = ebc.imb.se,
                                      EBC_imb_abs = NA,
                                      EBC_imb_abs_se = NA,
                                      EBC_imb_cum = NA,
                                      EBC_imb_cum_ciL = NA,
                                      EBC_imb_cum_ciU = NA,
                                      EBC_ratio = ebc.ratio,
                                      datacov_a = datacov_a,
                                      datacov_p = datacov_p
                                    )
        )
      } else {
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    data.frame(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "mensile",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var
                                    )
        )
      }
      
    } # regression loop over TE END
    
  } # regression loop over AE END
  
  cat('done.\n')  # scala stagionale
  
  
  # ---- Scala annuale  ----
  # EBC stats calcolate a scala annuale, VALORI AGGREGATI
  
  cat('Scala annuale (gap-filled):\n')
  
  # Converte timestamp in POSIXct
  # TS.cur <- c()
  # TS.cur <- EBC.base.data[[iebc]]$'TIMESTAMP_START'
  # TS.cur <- lubridate::ymd_hm(as.character(TS.cur))
  
  # Aggregazione annuale di AE
  # Format timestamps as year
  TS.yea <- format(TS.cur, '%Y')
  AE_year <- data.frame()
  AE_year <- AE %>%
    group_by(year = TS.yea) %>%
    summarise(across(everything(), sum, na.rm = TRUE))
  AE_year <- as.data.frame(AE_year)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  AE_year <- AE_year %>% mutate_if(is.numeric, ~na_if(., 0))
  
  # Aggregazione annuale di TE
  TE_year <- data.frame()
  TE_year <- TE %>%
    group_by(year = TS.yea) %>%
    summarise(across(everything(), sum, na.rm = TRUE))
  TE_year <- as.data.frame(TE_year)
  
  # Drop 0s caused by NAs in the half-hours (the period at the beginning)
  # Replace only on all Numeric columns
  TE_year <- TE_year %>% mutate_if(is.numeric, ~na_if(., 0))
  
  
  for (ae_var in names(AE)) { # regression loop over AE
    
    for (te_var in names(TE)) { # regression loop over TE 
      
      dati <- data.frame(NULL)
      dati <- data.frame(
        x = AE_year[[ae_var]],
        y = TE_year[[te_var]]
      )
      dati <- na.omit(dati)
      
      if (nrow(dati) > 0 & !any(apply(dati, 2, sum, na.rm=T) == 0)) {
        
        # EBC computation ........................
        
        # 1. Regressions: RMA and OLS (free and 0 intercept))
        if (nrow(dati) > 2) {  # compute regressions only if more then 2 years are available
          res <- data.frame()
          res <- calc_reg(dati$x, dati$y) 
        } else {
          cat('\nSlopes at annual scale cannot be computed: less then 3 years available!\n')
          res <- suppressWarnings(calc_reg(1:10, 1:10))
          res[1, ] <- NA
        }
        
        # 2. EBC Imbalance
        ebc.imb <- NULL; ebc.imb.se <- NULL 
        if (imb.as.mean) {
          ebc.imb <- mean(dati$x - dati$y, na.rm=T) # mean
        } else {
          ebc.imb <- median(dati$x - dati$y, na.rm=T) # median
        }  
        ebc.imb.se <- sd(dati$x - dati$y, na.rm = TRUE) / sqrt(length(na.omit(dati$x - dati$y)))
        
        # 3. EBC ratio: TE/AE
        ebc.ratio <- NULL
        ebc.ratio <- sum(dati$y, na.rm=TRUE) / sum(dati$x, na.rm=T)
        
        # Data coverage
        datacov_a <- nrow(dati)
        datacov_p <- round(nrow(dati)/nrow(AE_year)*100,2)
        
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    cbind(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "annuale",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var,
                                      res,
                                      EBC_imb = ebc.imb,
                                      EBC_imb_se = ebc.imb.se,
                                      EBC_imb_abs = NA,
                                      EBC_imb_abs_se = NA,
                                      EBC_imb_cum = NA,
                                      EBC_imb_cum_ciL = NA,
                                      EBC_imb_cum_ciU = NA,
                                      EBC_ratio = ebc.ratio,
                                      datacov_a = datacov_a,
                                      datacov_p = datacov_p
                                    )
        )
      } else {
        EBC_results_GF <- bind_rows(EBC_results_GF,
                                    data.frame(
                                      site = cur.site,
                                      network = cur.ntwrk,
                                      igbp = cur.igbp,
                                      t_scale = "mensile",
                                      fattore = NA,
                                      livello = NA,
                                      AE_var = ae_var,
                                      TE_var = te_var
                                    )
        )
      }
      
    } # regression loop over TE END
    
  } # regression loop over AE END
  
  cat('done.\n') # scala annuale
  
  
  # Execution time
  tx <- Sys.time()
  ellapsed.time <- difftime(tx, t0, units="hours")
  cat(paste('\n......', cur.site, 'EBC regressions computed in:', round(ellapsed.time, 3), 'hours ......\n\n'))
  
}    # loop over stations END


# Save it!
# summary(EBC_results_GF)
fwrite(EBC_results_GF, paste0('OUTPUT/EBC_results_gapfilled_v', vernum, '.csv'), sep=',', na = NA, quote = F, row.names = F)


# END ...

