#########################################################
# Net Ecosystem Exchange from Ameriflux flux tower data           #
#########################################################

#########################################################
# Santarem, Brazil                                       #
#########################################################


loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/BR-Sa1/BRSa12002_L4_h.txt"); # AMERIFLUX br_S1 Santarem Brazil
test=download.file(loc,"temp.dat");  
sa_2002=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/BR-Sa1/BRSa12003_L4_h.txt"); # AMERIFLUX br_S1 Santarem Brazil
test=download.file(loc,"temp.dat");  
sa_2003=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/BR-Sa1/BRSa12004_L4_h.txt"); # AMERIFLUX br_S1 Santarem Brazil
test=download.file(loc,"temp.dat");  
sa_2004=read.csv("temp.dat", header=TRUE);




santarem <- rbind(sa_2002,sa_2003,sa_2004); # merged dataframes


temp <- (santarem[,] == -9999.00); #undefined
santarem[temp]=NA;

N_obs <- length(santarem$DoY);
dt <- (2005-2002)/N_obs
t_obs <- seq(1992,2007-dt,by=dt);

NEE.santarem.ts <- ts(santarem$NEE_st_fMDS,start=2002,deltat=dt);

#night data

temp <- (santarem$Rg_f > 10 | is.na(santarem$Rg_f)); # day when Rg > 10W/m2
NEE.santarem.night <- santarem$NEE_st_fMDS;
NEE.santarem.night[temp]=NA;

NEE.santarem.night.ts <- ts(NEE.santarem.night,start=2002,deltat=dt);


# cumulative NEE time-series, treating missing data as zero NEE

temp <- is.na(santarem$NEE_st_fMDS);
santarem_cumul=santarem$NEE_st_fMDS;
santarem_cumul[temp] = 0;
santarem_cumul = cumsum(santarem_cumul)

NEE.santarem.cumulative.ts <- ts(santarem_cumul,start=2002,deltat=dt) #cumulative NEE time-series



#######################################
# Harvard forest deciduous broadeaf   #
#######################################

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11992_L4_h.txt"); # AMERIFLIX US_Ha1 Harvard Forest
test=download.file(loc,"temp.dat");  
ha_1992=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11993_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1993=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11994_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1994=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11995_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1995=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11996_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1996=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11997_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1997=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11998_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1998=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa11999_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_1999=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12000_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2000=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12001_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2001=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12002_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2002 =read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12003_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2003=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12004_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2004=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12005_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2005=read.csv("temp.dat", header=TRUE);

loc=file.path("ftp://cdiac.ornl.gov/pub/ameriflux/data/Level4/Sites_ByID/US-Ha1/USHa12006_L4_h.txt"); # AMERIFLIX US_Ha1 Santarem Brazil
test=download.file(loc,"temp.dat");  
ha_2006=read.csv("temp.dat", header=TRUE);


names(ha_1992);

# merged dataframes

harvard <- rbind(ha_1992,ha_1993,ha_1994,ha_1995,ha_1996,ha_1997,ha_1998,ha_1999,ha_2000,ha_2001,ha_2002,ha_2003,ha_2004,ha_2005,ha_2006);

#NAs

temp <- (harvard[,] == -9999.00);
harvard[temp]=NA;

#save data

#write.table(harvard,"e:\\redEdge\Biospherica\EddyCovariance\ha1Flux.dat",row.names=False)

#t_observations & NEE timeseries

N_obs <- length(harvard$DoY);
dt <- (2007-1992)/N_obs
t_obs <- seq(1992,2007-dt,by=dt);

NEE.harvard.ts <- ts(harvard$NEE_st_fMDS,start=1992,deltat=dt)

#night data harvard

temp <- (harvard$Rg_f > 10 | is.na(harvard$Rg_f)); # sunlight
NEE.harvard.night <- harvard$NEE_st_fMDS;
NEE.harvard.night[temp]=NA;

NEE.harvard.night.ts <- ts(NEE.harvard.night,start=1992,deltat=dt); #harvard NEE time-series


scale = 30*60*10^(-6)/1000*12  # converts to KgC/m2

temp <- is.na(harvard$NEE_st_fMDS);
harvard_cumul=harvard$NEE_st_fMDS;
harvard_cumul[temp] = 0;
harvard_cumul = cumsum(harvard_cumul);

NEE.harvard.cumulative.ts <- ts(harvard_cumul,start=1992,deltat=dt);


# NEE for both ecosystems with Night/Day distinction

par(mfrow=c(2,1))

plot( 12*NEE.harvard.ts/1000,xlab="", lwd=1,col=3,font.lab=2,font.axis=2,ylim=c(-0.6,0.5));
lines(12*NEE.harvard.night.ts/1000,lwd=2,col=1)
title(main="Harvard Forest CO2 Flux 1993-2006",cex.lab=2,cex.main=2)
title(ylab="mg C/m2/s",cex.lab=1)
legend("bottomleft", legend=c("Night","Day"), col=c(1,3),lwd=2:1,bty="n");
plot( 12*NEE.santarem.ts/1000, pch=16,xlab="", lwd=2,col=2,font.lab=2,font.axis=2,ylim=c(-0.6,0.5));
lines(12*NEE.santarem.night.ts/1000,lwd=2,col=1)
title(main="Tapajos Forest CO2 Flux 2002-2005",cex.lab=2,cex.main=2)
legend("bottomleft", legend=c("Night","Day"), col=c(1,2),lwd=2:1,bty="n");
text(2002.5,0.45, "data source http://public.ornl.gov/ameriflux/");
title(ylab="mg C/m2/s",cex.lab=1)
par(mfrow=c(1,1))

# plot cumulative NEE

plot( NEE.harvard.cumulative.ts*scale, pch=16,xlab="", lwd=2,col=3,font.lab=2,font.axis=2,ylim=c(-4.5,2),ylab="KgC/m2");
text(1996,-3, "data source http://public.ornl.gov/ameriflux/");
legend("bottomleft", legend=c("Harvard Forest 1993-2006","km 67 Santarem 2002-2005"), col=c(3,2),lwd=2:2,bty="n");
lines( NEE.santarem.cumulative.ts*scale, pch=16,xlab="",lwd=2,col=2,font.lab=2,font.axis=2);
title(ylab="KgC/m2",main="Cumulative NEE",cex.lab=1,cex.main=2)