########################################################
# calculate IOM in t/ha
# C input
# total carbon stock in t/ha
#####################################################

fIOM.Falloon.RothC =function(SOC, par1=-1.31, par2=1.139)
{
  
  # IOM=10^(par1+par2*log10(c))
  IOM=0.049*SOC^(1.139) 
  IOM
}
#################################################################################
# fget_equilibrium_fractions.RothC_input 
# brief: quantifies pool distribution and C input for RothC at equilibrium
#
#Input
# xi= scalar representing an averaged modifying factor
# C.tot = initial C stock (and C stock in equilibrium)
# clay = clay content
# fractI = vector of Cinput fractions that enter the DPM, RPM, HUM   
#          with a DR of 1.44 fractI becomes [1] 0.5901639 0.4098361 0.0000000
#          by fractI=c((DR)/(DR+1),1-(DR)/(DR+1),0)
#Output
# list with pools at equilibrium and C input at equilibrium
################################################################################
fget_equilibrium_fractions.RothC_input=function(xi=xi.frame,C.tot,clay, fractI)
{   
  rmf=xi
  IOM= fIOM.Falloon.RothC(SOC = C.tot)
  C.active=C.tot-IOM 
  
  ########################################################################
  #The analytical solution of RothC
  ########################################################################
  
  ########################################################################
  # Parameter
  ########################################################################
  fract.rooted.to.bio = 0.46
  fract.rooted.to.hum = 0.54
  ks = c(k.DPM = 10, k.RPM = 0.3, k.BIO = 0.66, k.HUM = 0.02, 
         k.IOM = 0)
  ks=as.numeric(ks)
  k.dpm=ks[1]
  k.rpm=ks[2]
  k.bio=ks[3]
  k.hum=ks[4]
  ########################################################################
  # the carbon use efficiency
  ########################################################################
  cue=  1/(1+ 1.67 * (1.85 + 1.6 * exp(-0.0786 * clay)))
  
  ########################################################################
  # All the coefficients alpha.1 und alpha.2
  ########################################################################
  alpha.1=cue*fract.rooted.to.bio
  alpha.2=cue*fract.rooted.to.hum
  
  ########################################################################
  # All the coefficients a.1.1, a.1.2, a.2.1, a2.2
  ########################################################################
  a.1.1=k.bio*rmf*(alpha.1-1)
  a.1.2=alpha.1*k.hum*rmf
  a.2.1=alpha.2*k.bio*rmf
  a.2.2=k.hum*rmf*(alpha.2-1)
  
  #########################################################################
  #########################################################################
  # The Eigenvalues lambda 1 and lambda 2
  #########################################################################
  lambda.1= (a.1.1+a.2.2)/2-sqrt(((a.1.1+a.2.2)/2)*((a.1.1+a.2.2)/2)+a.1.2*a.2.1-a.1.1*a.2.2)
  lambda.2= (a.1.1+a.2.2)/2+sqrt(((a.1.1+a.2.2)/2)*((a.1.1+a.2.2)/2)+a.1.2*a.2.1-a.1.1*a.2.2)
  #########################################################################
  # The c.0.1; c.0.2; c.0.3 values
  #########################################################################
  c.0.1= (alpha.2 * a.1.2 - alpha.1 * a.2.2)/(a.1.1*a.2.2-a.1.2*a.2.1)
  c.0.2= (alpha.2 * a.1.2 - alpha.1 * a.2.2)/(a.1.1*a.2.2-a.1.2*a.2.1)
  c.0.3= (a.1.2)/(a.1.1*a.2.2-a.1.2*a.2.1)
  
  ######################################################################################################
  # BIO pool quantification
  ######################################################################################################
  u.bio.dpm=(c.0.2) #65
  u.bio.rpm=(c.0.1) #66
  u.bio.hum=(c.0.3) #67
  
  
  ######################################################################################################
  # HUM pool quantification ( is all C.78)
  ######################################################################################################
  u.hum.dpm= 1/a.1.2*((-c.0.2*a.1.1-alpha.1))
  u.hum.rpm= 1/a.1.2*(-c.0.2*a.1.1-alpha.1)
  u.hum.hum= 1/a.1.2*(-c.0.3*a.1.1)
  
  
  ######################################################################################################
  # DPM C ( is all C.79)
  ######################################################################################################
  u.dpm.dpm=1/k.dpm/rmf 
  
  #C.dpm=i.dpm * u.dpm.dpm + C0 * s.dpm
  
  
  ######################################################################################################
  # RPM C ( is all C.80)
  ######################################################################################################
  u.rpm.rpm=1/k.rpm/rmf
  
  #C.rpm=i.rpm * u.rpm.rpm + C0 *s.rpm
  
  ######################################################################################################
  # Total C ( is all C.78)
  ######################################################################################################
  u.dpm=u.dpm.dpm+u.bio.dpm+u.hum.dpm
  u.rpm=u.rpm.rpm+u.bio.rpm+u.hum.rpm
  u.hum=u.bio.hum+u.hum.hum
  
  denominator= fractI[1]*u.dpm+fractI[2]*u.rpm+fractI[3]*u.hum
  
  fract.dpm= fractI[1]*u.dpm.dpm/denominator
  fract.rpm= fractI[2]*u.rpm.rpm/denominator
  fract.bio= (fractI[1]*u.bio.dpm+fractI[2]*u.bio.rpm+fractI[3]*u.bio.hum)/denominator
  fract.hum= (fractI[1]*u.hum.dpm+fractI[2]*u.hum.rpm+fractI[3]*u.hum.hum)/denominator   
  
  fract.all=c(fract.dpm,fract.rpm,fract.bio,fract.hum)
  
  ###################################################
  # so unfortunately we have the IOM
  ###################################################
  fract.all_stock=(fract.all*C.active)
  fract.all=fract.all_stock/C.tot
  fract.all=append(fract.all,IOM/C.tot)
  pools=fract.all*C.tot
  Cin=(C.tot-pools[5])/denominator
  list(pools,Cin)
}