% Woody vegetation carbon fluxes

% NPP: calculate inputs from NPP to living pools
leafinput           = NPP(:,month) ./ 3;
woodinput           = NPP(:,month) ./ 3;
frootinput          = NPP(:,month) ./ 3;

% NPP: transfer NPP into living biomass pools
leafpool          = leafpool      + leafinput;
abovewoodpool     = abovewoodpool + woodinput * aboveWoodFraction;
belowwoodpool     = belowwoodpool + woodinput * (1 - aboveWoodFraction);
frootpool         = frootpool     + frootinput;

% HERBIVORY
herbivory           = trees_herbivory .* herb_seasonality(:,month);         % yearly herbivory * seasonality scalar
io                  = herbivory > leafpool;                                 % check that herbivory does not exceed available leaf
herbivory(io)       = leafpool(io);                                         % in case herbivory exceeds leaf, lower herbivory
leafpool            = leafpool - herbivory;                                 % deduct herbivory from the leafpool
carbonout_leaf      = herbivory .* (1 - herbivoreEff);                      % part of the consumed leaf will be returned as litter
herbivory           = herbivory  - herbivory .* (1 - herbivoreEff);         % part of the consumed leaf for maintenance

surfstrpool         = surfstrpool   + carbonout_leaf .*(1 - metabfract);    
surfmetpool         = surfmetpool   + carbonout_leaf .* metabfract;    

% DECAY of biomass and litter, each of the following equations have the following basic form: 
% carbon pool size * rate constant * abiotic effect. Some may have more terms but all are first order
carbonout_leaf      = leafpool      .* annK_leaf    .* litterscalar(:,month);
carbonout_abovewood = abovewoodpool .* K_wood;
carbonout_belowwood = belowwoodpool .* K_wood;
carbonout_froot     = frootpool     .* annK_froot   .* rootlitscalar(:,month);
carbonout_cwd       = cwdpool       .* K_cwd        .* a(:,month);
carbonout_surfmet   = surfmetpool   .* K_surfmet    .* a(:,month);
carbonout_surfstr   = surfstrpool   .* K_surfstr    .* a(:,month)  .* lignineffect;
carbonout_soilmet   = soilmetpool   .* K_soilmet    .* a(:,month);
carbonout_soilstr   = soilstrpool   .* K_soilstr    .* a(:,month)  .* lignineffect;
carbonout_surfmic   = surfmicpool   .* K_surfmic    .* a(:,month);
carbonout_soilmic   = soilmicpool   .* K_soilmic    .* a(:,month)  .* soilmicDecayFactor;
carbonout_slow      = slowpool      .* K_slow       .* a(:,month);
carbonout_armored   = armoredpool   .* K_armored    .* a(:,month);

% determine inputs into structural and metabolic pools from decaying living
% pools
surfstrpool       = surfstrpool   + (carbonout_leaf  + carbonout_cwd)       .*(1 - metabfract);    
soilstrpool       = soilstrpool   + (carbonout_froot + carbonout_belowwood) .*(1 - metabfract);
surfmetpool       = surfmetpool   + (carbonout_leaf  + carbonout_cwd)       .* metabfract;    
soilmetpool       = soilmetpool   + (carbonout_froot + carbonout_belowwood) .* metabfract;    
cwdpool           = cwdpool       + carbonout_abovewood;

leafpool          = leafpool      - carbonout_leaf;
abovewoodpool     = abovewoodpool - carbonout_abovewood;
belowwoodpool     = belowwoodpool - carbonout_belowwood;
frootpool         = frootpool     - carbonout_froot;
cwdpool           = cwdpool       - carbonout_cwd;
surfstrpool       = surfstrpool   - carbonout_surfstr;

% empty respirationpools in beginning of month
resppool_surfstr    = zeros(numberOfGridCells,1);
resppool_surfmet    = zeros(numberOfGridCells,1);
resppool_surfmic    = zeros(numberOfGridCells,1);
resppool_soilstr    = zeros(numberOfGridCells,1);
resppool_soilmet    = zeros(numberOfGridCells,1);
resppool_soilmic    = zeros(numberOfGridCells,1);
resppool_slow       = zeros(numberOfGridCells,1);
resppool_armored    = zeros(numberOfGridCells,1);

% respiratory fluxes from every pool (temp)
temp              = (carbonout_surfstr    .* structuralLignin)            .* eff_surfstr2slow;
slowpool          = slowpool              + temp;
resppool_surfstr  = resppool_surfstr      + (temp ./ eff_surfstr2slow)    .* (1 - eff_surfstr2slow);

temp              = (carbonout_surfstr    .* (1 - structuralLignin))      .* eff_surfstr2surfmic;
surfmicpool       = surfmicpool           + temp;
resppool_surfstr  = resppool_surfstr      + (temp ./ eff_surfstr2surfmic) .* (1 - eff_surfstr2surfmic);
soilstrpool       = soilstrpool           - carbonout_soilstr;

temp              = carbonout_soilstr     .* structuralLignin             .* eff_soilstr2slow;
slowpool          = slowpool               + temp;
resppool_soilstr  = resppool_soilstr       + (temp ./ eff_soilstr2slow)    .* (1 - eff_soilstr2slow);

temp              = carbonout_soilstr     .* (1 - structuralLignin)       .* eff_soilstr2soilmic;
soilmicpool       = soilmicpool            + temp;
resppool_soilstr  = resppool_soilstr       + (temp ./ eff_soilstr2soilmic) .* (1 - eff_soilstr2soilmic);

temp              = carbonout_surfmet      .* eff_surfmet2surfmic;
surfmetpool       = surfmetpool            - carbonout_surfmet;
surfmicpool       = surfmicpool            + temp;
resppool_surfmet  = (temp ./ eff_surfmet2surfmic) .* (1 - eff_surfmet2surfmic);

temp              = carbonout_soilmet      .* eff_soilmet2soilmic;
soilmetpool       = soilmetpool            - carbonout_soilmet;
soilmicpool       = soilmicpool            + temp;
resppool_soilmet  = (temp ./ eff_soilmet2soilmic) .* (1 - eff_soilmet2soilmic);

temp              = carbonout_surfmic      .* eff_surfmic2slow;
surfmicpool       = surfmicpool            - carbonout_surfmic;
slowpool          = slowpool               + temp;
resppool_surfmic  = (temp ./ eff_surfmic2slow)    .* (1 - eff_surfmic2slow);

resppool_soilmic  = eff_soilmic2slow       .* carbonout_soilmic;
soilmicpool       = soilmicpool            - carbonout_soilmic;

temp              = carbonout_soilmic      .* (0.003 + (0.032 .* clay));
armoredpool       = armoredpool            + temp;

temp              = carbonout_soilmic      - temp   - resppool_soilmic;
slowpool          = slowpool              + temp;

resppool_slow     = carbonout_slow        .* (1 - eff_slow2soilmic);
slowpool          = slowpool               - carbonout_slow;

temp              = carbonout_slow        .* eff_slow2soilmic .* decayClayFactor;
armoredpool       = armoredpool          + temp;

temp              = carbonout_slow        - resppool_slow     - temp;
soilmicpool       = soilmicpool          + temp;

temp              = carbonout_armored     .* eff_armored2soilmic;
armoredpool       = armoredpool          - carbonout_armored;
soilmicpool       = soilmicpool          + temp;

resppool_armored  = (temp ./ eff_armored2soilmic) .* (1 - eff_armored2soilmic);

% FIRES consume part of the pools depending on burned fraction (BF),
% combustion completeness (CC), and tree mortality rate 
combusted_leaf          = leafpool        .* BF(:,month)  .* ccLeaf(:,month)      .* mortality_tree;
combusted_abovewood     = abovewoodpool   .* BF(:,month)  .* ccWood(:,month)      .* mortality_tree;
combusted_cwd           = cwdpool         .* BF(:,month)  .* ccCwd(:,month);
combusted_surfstr       = surfstrpool     .* BF(:,month)  .* ccFineLitter(:,month);
combusted_surfmet       = surfmetpool     .* BF(:,month)  .* ccFineLitter(:,month);
combusted_surfmic       = surfmicpool     .* BF(:,month)  .* ccFineLitter(:,month);

% FIRE: the non-combusted parts
nonCombusted_leaf       = leafpool        .* BF(:,month)  .* (ones(numberOfGridCells,1) - ccLeaf(:,month)) .* mortality_tree;
nonCombusted_abovewood  = abovewoodpool   .* BF(:,month)  .* (ones(numberOfGridCells,1) - ccWood(:,month)) .* mortality_tree;
nonCombusted_belowwood  = belowwoodpool   .* BF(:,month)  .* mortality_tree;
nonCombusted_froot      = frootpool       .* BF(:,month)  .* mortality_tree;

% FIRE flux from not combusted parts to other pools
surfstrpool           = surfstrpool +  nonCombusted_leaf .* (1 - metabfract);
surfmetpool           = surfmetpool +  nonCombusted_leaf .* metabfract;   
soilstrpool           = soilstrpool + (nonCombusted_froot + nonCombusted_belowwood) .* (1 - metabfract);
soilmetpool           = soilmetpool + (nonCombusted_froot + nonCombusted_belowwood) .* metabfract;   
cwdpool               = cwdpool     +  nonCombusted_abovewood;

% FIRE 
leafpool              = leafpool        - combusted_leaf        - nonCombusted_leaf;
abovewoodpool         = abovewoodpool   - combusted_abovewood   - nonCombusted_abovewood;
belowwoodpool         = belowwoodpool                           - nonCombusted_belowwood;
frootpool             = frootpool                               - nonCombusted_froot;
cwdpool               = cwdpool         - combusted_cwd;
surfstrpool           = surfstrpool     - combusted_surfstr;
surfmetpool           = surfmetpool     - combusted_surfmet;
surfmicpool           = surfmicpool     - combusted_surfmic;

% FUELWOOD collection
if numberOfAgeClasses == 1
	fuelwoodout         = fuelwooddemand + fuelshortage;                        % fuel wood demand 
	io                  = fuelwoodout > cwdpool;                                % in case demand exceeds availability
	fuelwoodout(io)     = cwdpool(io);                                          % demand becomes availability
	fuelshortage(io)    = fuelshortage(io) - fuelwoodout(io) + cwdpool(io);     % and shortage increases
	io                  = fuelwoodout < cwdpool;                                % in case availability exceeds demand
	fuelshortage(io)    = fuelshortage(io) - cwdpool(io) + fuelwoodout(io);     % shortage decreases
	io                  = fuelshortage < 0;
	fuelshortage(io)    = 0;
	cwdpool             = cwdpool - fuelwoodout;                                % fuel wood is taken out of the cwd pool
else
    fuelwoodout         = fuelwooddemand;
    io                  = fuelwooddemand > cwdpool;
    fuelwoodout(io)     = cwdpool(io);
    cwdpool             = cwdpool - fuelwoodout;
end
    
% calculate fluxes
if numberOfAgeClasses==1
	wresp          = resppool_surfstr + resppool_surfmet + resppool_surfmic + resppool_soilstr + ...
                                resppool_soilmet + resppool_soilmic + resppool_slow + resppool_armored;
	wcomb          = combusted_leaf + combusted_abovewood + combusted_cwd + combusted_surfstr + ...
                                combusted_surfmet + combusted_surfmic;
	wherb          = herbivory;
	wbiof          = fuelwoodout;
else
    if ageClass == 1
        wresp = zeros(numberOfGridCells,1);
        wcomb = zeros(numberOfGridCells,1);
        wherb = zeros(numberOfGridCells,1); 
        wbiof = zeros(numberOfGridCells,1);
    end
	wresp          = wresp + (resppool_surfstr + resppool_surfmet + resppool_surfmic + resppool_soilstr + resppool_soilmet + ...
                                resppool_soilmic + resppool_slow + resppool_armored) ./ numberOfAgeClasses;
	wcomb          = wcomb + (combusted_leaf + combusted_abovewood + combusted_cwd + combusted_surfstr + combusted_surfmet + ...
                                combusted_surfmic) ./ numberOfAgeClasses;
	wherb          = wherb + herbivory ./ numberOfAgeClasses;
	wbiof          = wbiof + fuelwoodout ./ numberOfAgeClasses;
end    
% *************************************************************************
