source('initvars.R')
source('AllocateMo.R')
source('AllocateYr.R')
source('AtmEnviron.R')
source('CNTrans.R')
source('Decomp.R')
source('Leach.R')
source('Phenology.R')
source('Photosyn.R')
source('SoilResp.R')
source('storeoutput.R')
source('Waterbal.R')
source('YearInit.R')

spinup_pnetcn_continue <- function(clim,site,veg,share,reps)
{
  yearnum=1+max(clim$year)-min(clim$year)
  climtemp=clim
  climtemp$year=climtemp$year-yearnum*(reps-1)
  share=initvars(climtemp,veg)
  share$FolMassMax=NA
  share$FolMassMin=NA
  share$FolNCon=NA
  share$year=climtemp$year[1]
  share$doy=climtemp$doy[1]
  share$time.POSIXt=strptime(paste(share$year,share$doy,sep="-"),format="%Y-%j")
  outyearly=share
  outmonthly=share
  for (j in 1:reps)
  {
    output=pnetcn_continue(climtemp,site,veg,share)
    print(paste("length(outmonthly)=",length(outmonthly)))
    print(paste("length(outyearly)=",length(outyearly)))
    print(paste("length(output[1:(length(output)/2)])=",length(output[1:(length(output)/2)])))
    newoutmonthly=as.data.frame(output[1:(length(output)/2)])
    outmonthly=rbind(as.data.frame(outmonthly),newoutmonthly[2:length(newoutmonthly[,1]),])
    outyearly=rbind(as.data.frame(outyearly),as.data.frame(output[(1+length(output)/2):length(output)]))
    share=outyearly[length(outyearly[,1]),]
    climtemp$year=climtemp$year+yearnum
  }
  return(c(outmonthly,outyearly));
  #return(output);
}

spinup_clim_calcCO2 <- function(clim,reps){
  print("num1")
  yearnum=1+max(clim$year)-min(clim$year)
  print(head(clim$year))
  climnew=clim
  climnew$year=clim$year-yearnum*(reps-1)
  print(paste("yearnum=",yearnum,"reps=",reps))
  print(head(climnew$year))
  climtemp=climnew
  print("num2")
  SaveCO2=climtemp$CO2
  print(paste("head(SaveCO2)",head(SaveCO2)))
  SaveCO2[ is.na(clim$CO2)]=CO2Fun(as.numeric(clim$year[ is.na(clim$CO2)]))
  print(paste("head(SaveCO2)",head(SaveCO2)))
  print("num3")
  #SaveCO2=climtemp$CO2
  print(paste("head(climnew$year)=",head(climnew$year)))
  climnew$CO2=CO2Fun(climnew$year) #need to check which of these is which
  print(paste("head(climnew$CO2)=",head(climnew$CO2)))
  #climtemp$CO2=CO2Fun(climtemp$year) #need to check which of these is which
  for (i in 1:(reps-1)){
    climtemp$year=climtemp$year+yearnum
    if(i==(reps-1)){
      print("num4")
      climtemp$CO2=SaveCO2
      print("num5")
    }else{
      climtemp$CO2=CO2Fun(climtemp$year)
    }
    climnew=rbind(climnew,climtemp)
  }
  return(climnew)
}

CO2Fun=function(yearVector){
  #Based on: clim->CO2[i] = 280 +( (clim->year[i]-1800)*0.0188)^3.35
  #from Scott Ollinger, 2002,GCB,8,545-562
  #see also: "CO2_FromZaixing.xlsx"
  print(length(is.na(yearVector)))
  print(head(yearVector,40))
  CO2=rep(280,length(yearVector))
  CO2[yearVector>1800]=280 + ((yearVector[yearVector>1800]-1800)*0.0188)^3.35
  return(CO2)
}

spinup_clim <- function(clim,reps){
  yearnum=1+max(clim$year)-min(clim$year)
  climnew=clim
  climnew$year=clim$year-yearnum*(reps-1)
  climtemp=climnew
  for (i in 1:(reps-1)){
    climtemp$year=climtemp$year+yearnum
    climnew=rbind(climnew,climtemp)
  }
  return(climnew)
}

spinup_clim2 <- function(clim1,reps1,clim2,reps2){
  #Assumes clim1 and clim2 end on the same year
  yearnum1=1+max(clim1$year)-min(clim1$year)
  yearnum2=1+max(clim2$year)-min(clim2$year)
  climnew1=clim1
  #climnew1$year=clim1$year-yearnum1*(reps1-1)-yearnum2*(reps2-1)
  climnew1$year=clim1$year-yearnum1*(reps1-1)-yearnum2*(reps2)
  climtemp1=climnew1
  for (i in 1:(reps1-1)){
    climtemp1$year=climtemp1$year+yearnum1
    climnew1=rbind(climnew1,climtemp1)
  }
  climnew2=clim2
  climnew2$year=clim2$year-yearnum2*(reps2-1)
  climtemp2=climnew2
  for (i in 1:(reps2-1)){
    climtemp2$year=climtemp2$year+yearnum2
    climnew2=rbind(climnew2,climtemp2)
  }
  climnew=rbind(climnew1,climnew2)
  return(climnew)
}

spinup_pnetcn <- function(clim,site,veg,reps,toDate=TRUE)
{
  yearnum=1+max(clim$year)-min(clim$year)
  climtemp=clim
  if(toDate==TRUE)
  {
    climtemp$year=climtemp$year-yearnum*(reps-1)
  }else{
    climtemp$year=climtemp$year-yearnum*(reps)
  }
  share=initvars(climtemp,veg)
  share$FolMassMax=NA
  share$FolMassMin=NA
  share$FolNCon=NA
  share$year=climtemp$year[1]
  share$doy=climtemp$doy[1]
  share$TransMo=NA
  share$time.POSIXt=strptime(paste(share$year,share$doy,sep="-"),format="%Y-%j")
  outyearly=share
  outmonthly=share
  for (j in 1:reps)
  {
    output=pnetcn_continue(climtemp,site,veg,share)
    print(paste("length(outmonthly)=",length(outmonthly)))
    print(paste("length(outyearly)=",length(outyearly)))
    print(paste("length(output[1:(length(output)/2)])=",length(output[1:(length(output)/2)])))
    newoutmonthly=as.data.frame(output[1:(length(output)/2)])
    outmonthly=rbind(as.data.frame(outmonthly),newoutmonthly[2:length(newoutmonthly[,1]),])
    outyearly=rbind(as.data.frame(outyearly),as.data.frame(output[(1+length(output)/2):length(output)]))
    share=outyearly[length(outyearly[,1]),]
    climtemp$year=climtemp$year+yearnum
  }
  return(c(outmonthly,outyearly));
  #return(output);
}

pnetcn_continue <-function(climIn,siteIn,vegIn,shareIn) {  
  #
  # [out]<-pnetcn(climIn,siteIn,vegIn)
  #
  
  ## Initialize
  #clear global
  #
  ## Vegetation and site constants
  #global veg     # Vegetation parameters (input)
  #global site    # Site specif (ic parameters (input)) {
  #global clim    # Climatological time series (input)
  #global share   # Internal variables shared between functions
  #global out     # Calculated variables int}#ended to be returned to the user
  #global rstep   # The current running timestep of the execution (day or month)
  #global ystep   # The yearstep (only used in constructing the output array)
  #global CN_Mode # 1 if ( running PnET-CN) {
  
  # Assign input structures to those used within the model
  veg <- vegIn  
  site <- siteIn  
  clim <- climIn
  share <- shareIn
  
  ystep=1
  rstep=1
  # Initialize variables
  share$year=clim$year[1]
  share$doy=clim$doy[1]
  share$time.POSIXt=strptime(paste(share$year,share$doy,sep="-"),format="%Y-%j")
  share=YearInit(share)
  outyearly=share
  outmonthly=share
  #outyearly$NetPsnMoNoO3=NA
  ##outyearly$SnowPack=NA
  #outmonthly$NetPsnMoNoO3=NA
  ##outmonthly$SnowPack=NA #APPARENTLY this modifies share... don't want
  print(paste("share$GDDTot=",share$GDDTot))
  clim$dayspan=clim$doy #just to get the length correct, real def. below
  #clim$dayspan=30
  for (  i in 1 : length(clim$doy) ) {
    clim$dayspan[i]<-30
  }#end
  
  CN_Mode<-1
  
  ## Main run loop
  # Loop through all climate records
  print(paste("rstep(1)=",rstep))
  print(paste("clim$doy[rstep](1)=",clim$doy[rstep]))
  for (  rstep in 1 : length(clim$doy) ) {
    print(paste("rstep(2)=",rstep))
    print(paste("clim$doy[rstep](2)=",clim$doy[rstep]))
    print(paste("(rstep!=1) =",rstep!=1))
  
    # End-of-year activity
    if ( (rstep != 1) && (clim$doy[rstep] < clim$doy[rstep-1])) {
      share$FolMassMax=veg$FolMassMax #to be able to change
      share$FolMassMin=veg$FolMassMin #to be able to change
      share$FolNCon=veg$FolNCon #to be able to change
      share=AllocateYr(share,clim,site,veg,rstep,1)
      print("I am here?")
      print(paste("veg$FolMassMax=", veg$FolMassMax))
      print(paste("share$FolMassMax=", share$FolMassMax))
      veg$FolMassMax=share$FolMassMax #changing ************
      veg$FolMassMin=share$FolMassMin #changing
      veg$FolNCon=share$FolNCon #changing
      #out=storeoutput(1)
      outyearly=rbind(outyearly,share)
      ystep=ystep+1
      share=YearInit(share)
      print(paste("rstep(3)=",rstep))
      print(paste("clim$doy[rstep](3)=",clim$doy[rstep]))
    }#end
  
    # Call subroutines
    print(paste("share$PlantN (1) (pnetcn, rstep=",rstep,") =",share$PlantN))
    share=AtmEnviron(share,clim,rstep)
    share=Phenology(share,clim,veg,rstep,1)
    share=Photosyn(share,clim,site,veg,rstep)
    share$SnowPack=site$SnowPack #so can be modified by Waterbal routine
    print(paste("share$SnowPack (before, rstep=",rstep,") =",share$SnowPack))
    print(paste("site$SnowPack (before, rstep=",rstep,") =",site$SnowPack))
    share=Waterbal(share,clim,site,veg,rstep)
    site$SnowPack=share$SnowPack #incorporates change to SnowPack
    print(paste("share$SnowPack (after, rstep=",rstep,") =",share$SnowPack))
    print(paste("site$SnowPack (after, rstep=",rstep,") =",site$SnowPack))
    #AllocateMo_NoTMult
    print(paste("share$PlantN (2) (pnetcn, rstep=",rstep,") =",share$PlantN))
    share=AllocateMo(share,clim,veg,rstep,CN_Mode)
    share=Phenology(share,clim,veg,rstep,2)
    print(paste("share$PlantN (2.5) (pnetcn, rstep=",rstep,") =",share$PlantN))
    share=CNTrans(share,clim,site,veg,rstep)
    print(paste("share$PlantN (3) (pnetcn, rstep=",rstep,") =",share$PlantN))
    share=Decomp(share,clim,veg,rstep)
    share=Leach(share)
    #out=storeoutput(outmonthly,share,rstep,0)
    #print(names(outmonthly))
    #print(names(share))
    #add dates
    share$year=clim$year[rstep]
    share$doy=clim$doy[rstep]
    share$time.POSIXt=strptime(paste(share$year,share$doy,sep="-"),format="%Y-%j")
    #print(names(outmonthly))
    #print(names(share))
    outmonthly=rbind(outmonthly,share)
  
  }#end
  
  # Calculate final year
  #AllocateYr
  share=AllocateYr(share,clim,site,veg,rstep,1)
  #storeoutput(1)
  outyearly=rbind(outyearly,share)
  return(c(outmonthly,outyearly));
}