## set necessary libraries
#require(fields)
#require(ncdf4)
#require(ncdf)
#require(M3)
library(stringr)

infile_txt <- Sys.getenv("INFILE")
outfile_csv <- Sys.getenv("OUTFILE")
nc_csv <- Sys.getenv("NCFILE")
refile_txt <- Sys.getenv("REFILE")
print(infile_txt)
print(outfile_csv)
print(refile_txt)

#source("func_for_CMAQ_files_with_ncdf_rgdal.r")

in.df <- data.frame(read.table(infile_txt,header=TRUE, skip=0, na.strings="NA", strip.white=TRUE))
re.df <- data.frame(read.table(refile_txt,header=TRUE, skip=0, na.strings="NA", strip.white=TRUE))

#Tracer site ssray  mon   day  year  fsrh  flrh fssrh   frh
#  date     NO2(ppm)   Soil          Amm_NO3        OMC
#  LAC      CM         Amm_SO4       Sea_Salt

# Add facility lines
for ( item in re.df[,1] ) {
  print(item)   # first item is F001
  print(re.df[,3][re.df[,1]==item])
  slist <- unlist(strsplit(as.vector(re.df[,3][re.df[,1]==item]), "\\,"))
  slist <- str_pad(slist, 3, pad = "0")
  slist <- paste("S",slist, sep="")
  # [1] "S002" "S003" "S004", take the first one as template
  indic <- in.df$Tracer == slist[1]
  newre.df <- in.df[indic,]
  newre.df$Tracer <- item

  # number of source length(slist) 
  ns <- length(slist)
  if ( ns > 1)      { 
  for ( i in 2:ns ) {
    indic <- in.df$Tracer == slist[i]
    temre.df <- in.df[indic,]
    newre.df$NO2   <- newre.df$NO2   + temre.df$NO2
    newre.df$Soil  <- newre.df$Soil  + temre.df$Soil
    newre.df$Amm_NO3 <- newre.df$Amm_NO3 + temre.df$Amm_NO3
    newre.df$OMC     <- newre.df$OMC   + temre.df$OMC 
    newre.df$LAC     <- newre.df$LAC   + temre.df$LAC 
    newre.df$CM      <- newre.df$CM    + temre.df$CM  
    newre.df$Amm_SO4 <- newre.df$Amm_SO4 + temre.df$Amm_SO4
    newre.df$Sea_Salt  <- newre.df$Sea_Salt + temre.df$Sea_Salt
  }  
  }
  in.df <- rbind(in.df, newre.df)
}


# Calculate large and small ( Adjusted concentration ug/m3 )
in.df$L_oc  <- ifelse (in.df$OMC<20.0,(in.df$OMC/20.0)*in.df$OMC, in.df$OMC)
in.df$S_oc  <- in.df$OMC - in.df$L_oc
in.df$L_so4 <- ifelse (in.df$Amm_SO4<20.0,(in.df$Amm_SO4/20.0)*in.df$Amm_SO4, in.df$Amm_SO4)
in.df$S_so4 <- in.df$Amm_SO4 - in.df$L_so4
in.df$L_no3 <- ifelse (in.df$Amm_NO3<20.0,(in.df$Amm_NO3/20.0)*in.df$Amm_NO3, in.df$Amm_NO3)
in.df$S_no3 <- in.df$Amm_NO3 - in.df$L_no3

# Light Extinction (Mm-1) 
in.df$E_oc  <- 2.8*in.df$S_oc + 6.1*in.df$L_oc
in.df$E_so4 <- 2.2*in.df$fsrh*in.df$S_so4 + 4.8*in.df$flrh*in.df$L_so4
in.df$E_no3 <- 2.4*in.df$fsrh*in.df$S_no3 + 5.1*in.df$flrh*in.df$L_no3
in.df$E_ec  <- 10*in.df$LAC
in.df$E_soil <- in.df$Soil
in.df$E_cm   <- 0.6*in.df$CM
in.df$E_no2  <- 0.33*in.df$NO2
in.df$E_salt <- 1.7*in.df$fssrh*in.df$Sea_Salt
in.df$T_bext <- in.df$E_oc+in.df$E_so4+in.df$E_no3+in.df$E_ec + in.df$E_soil + in.df$E_cm+in.df$E_no2
in.df$T_bext[in.df$Tracer=="Base_2016"] <- in.df$T_bext[in.df$Tracer=="Base_2016"]+in.df$E_salt[in.df$Tracer=="Base_2016"]+in.df$ss_rayleigh[in.df$Tracer=="Base_2016"]
in.df$dv <- 10*log(in.df$T_bext/10 )   # AN4/10 


# Light Extinction (Mm-1) by Scenario in the Absence of Source S**
# Ignore them right now
#siteNames <- unique(in.df$site)
#tracerID  <- unique(in.df$Tracer)
str(in.df)

# read Natural_Conditions_tabulated_fromMF.csv
#Class I Area,(1/Mm),dv
#BAND,14.00281,3.366726
nc.df <- data.frame(read.csv(nc_csv,header=TRUE, skip=0, na.strings="NA", strip.white=TRUE))
in.df$nc_ext <- 0.0
for ( item in nc.df[,1] ) {
  print(item)
  print(nc.df[,2][nc.df[,1]==item])
  in.df$nc_ext[grepl(item,in.df$site)] <- nc.df[,2][nc.df[,1]==item]
}
in.df$elta_dv <- 10*log((in.df$nc_ext+in.df$T_bext)/10)-10*log(in.df$nc_ext/10) 

# add percent: ratio of the tracer species to the base species at the same location
#siteNames <- unique(in.df$Tite)
tracerIDs <- unique(in.df$Tracer)
in.df$E_oc_p  <- 0.0
in.df$E_so4_p  <- 0.0
in.df$E_no3_p  <- 0.0
in.df$E_ec_p   <- 0.0
in.df$E_soil_p <- 0.0
in.df$E_cm_p   <- 0.0
in.df$E_no2_p  <- 0.0
in.df$T_bext_p  <- 0.0
in.df$dv_p  <- 0.0

indic <- in.df$Tracer == "Base_2016"
base.df <- in.df[indic,]
for ( item in tracerIDs ) {
  in.df$E_oc_p[in.df$Tracer==item]  <- (in.df$E_oc[in.df$Tracer==item]/base.df$E_oc)*100
  in.df$E_so4_p[in.df$Tracer==item] <- (in.df$E_so4[in.df$Tracer==item]/base.df$E_so4)*100
  in.df$E_no3_p[in.df$Tracer==item] <- (in.df$E_no3[in.df$Tracer==item]/base.df$E_no3)*100
  in.df$E_ec_p[in.df$Tracer==item]  <- (in.df$E_ec[in.df$Tracer==item]/base.df$E_ec)*100
  in.df$E_soil_p[in.df$Tracer==item] <- (in.df$E_soil[in.df$Tracer==item]/base.df$E_soil)*100
  in.df$E_cm_p[in.df$Tracer==item]   <- (in.df$E_cm[in.df$Tracer==item]/base.df$E_cm)*100
  in.df$E_no2_p[in.df$Tracer==item]  <- (in.df$E_no2[in.df$Tracer==item]/base.df$E_no2)*100
  in.df$T_bext_p[in.df$Tracer==item] <- (in.df$T_bext[in.df$Tracer==item]/base.df$T_bext)*100
  in.df$dv_p[in.df$Tracer==item]  <- (in.df$dv[in.df$Tracer==item]/base.df$dv)*100
}


write.table(in.df,outfile_csv,col.names=T,row.names=F, quote=FALSE,sep=",")
#write.table(format(in.df, digits=6),outfile_csv,col.names=T,row.names=F, quote=FALSE,sep=",")