% herbaceous vegetation carbon fluxes

% NPP: calculate inputs from NPP to living pools
leafinput           = NPP(:,month) .* 0.4;  % 40% of NPP goes aboveground, rest belowground
frootinput          = NPP(:,month) .* 0.6;

% NPP: transfer NPP into living biomass pools
hleafpool          = hleafpool      + leafinput;
hfrootpool         = hfrootpool     + frootinput;

% HERBIVORY
herbivory           = grass_herbivory .* herb_seasonality(:,month);         % yearly herbivory * seasonality scalar
io                  = herbivory > hleafpool;                                 % check that herbivory does not exceed available leaf
herbivory(io)       = hleafpool(io);                                         % in case herbivory exceeds leaf, lower herbivory
hleafpool           = hleafpool - 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

hsurfstrpool         = hsurfstrpool   + carbonout_leaf .*(1 - metabfract);    
hsurfmetpool         = hsurfmetpool   + 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      = hleafpool      .* K_hleaf      .* hlitterscalar(:,month);
carbonout_froot     = hfrootpool     .* K_hfroot     .* hlitterscalar(:,month);
carbonout_surfmet   = hsurfmetpool   .* K_surfmet    .* a(:,month);
carbonout_surfstr   = hsurfstrpool   .* K_surfstr    .* a(:,month)  .* lignineffect;
carbonout_soilmet   = hsoilmetpool   .* K_soilmet    .* a(:,month);
carbonout_soilstr   = hsoilstrpool   .* K_soilstr    .* a(:,month)  .* lignineffect;
carbonout_surfmic   = hsurfmicpool   .* K_surfmic    .* a(:,month);
carbonout_soilmic   = hsurfmicpool   .* K_soilmic    .* a(:,month)  .* soilmicDecayFactor;
carbonout_slow      = hslowpool      .* K_slow       .* a(:,month);
carbonout_armored   = harmoredpool   .* K_armored    .* a(:,month);

% determine inputs into structural and metabolic pools from decaying living
% pools
hsurfstrpool       = hsurfstrpool   + carbonout_leaf  .* (1 - metabfract);    
hsoilstrpool       = hsoilstrpool   + carbonout_froot .* (1 - metabfract);
hsurfmetpool       = hsurfmetpool   + carbonout_leaf  .* metabfract;    
hsoilmetpool       = hsoilmetpool   + carbonout_froot .* metabfract;    

hleafpool          = hleafpool      - carbonout_leaf;
hfrootpool         = hfrootpool     - carbonout_froot;
hsurfstrpool       = hsurfstrpool   - 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;
hslowpool          = hslowpool            + temp;
resppool_surfstr   = resppool_surfstr      + (temp ./ eff_surfstr2slow)    .* (1 - eff_surfstr2slow);

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

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

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

temp               = carbonout_surfmet     .* eff_surfmet2surfmic;
hsurfmetpool       = hsurfmetpool          - carbonout_surfmet;
hsurfmicpool       = hsurfmicpool          + temp;
resppool_surfmet   = (temp ./ eff_surfmet2surfmic) .* (1 - eff_surfmet2surfmic);

temp               = carbonout_soilmet     .* eff_soilmet2soilmic;
hsoilmetpool       = hsoilmetpool          - carbonout_soilmet;
hsurfmicpool       = hsurfmicpool           + temp;
resppool_soilmet   = (temp ./ eff_soilmet2soilmic) .* (1 - eff_soilmet2soilmic);

temp               = carbonout_surfmic     .* eff_surfmic2slow;
hsurfmicpool       = hsurfmicpool           - carbonout_surfmic;
hslowpool          = hslowpool              + temp;
resppool_surfmic   = (temp ./ eff_surfmic2slow)    .* (1 - eff_surfmic2slow);

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

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

temp               = carbonout_soilmic     - temp   - resppool_soilmic;
hslowpool          = hslowpool            + temp;

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

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

temp               = carbonout_slow        - resppool_slow     - temp;
hsurfmicpool       = hsurfmicpool          + temp;

temp               = carbonout_armored     .* eff_armored2soilmic;
harmoredpool       = harmoredpool          - carbonout_armored;
hsurfmicpool       = hsurfmicpool          + 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          = hleafpool        .* BF(:,month)  .* ccLeaf(:,month);
combusted_surfstr       = hsurfstrpool     .* BF(:,month)  .* ccFineLitter(:,month);
combusted_surfmet       = hsurfmetpool     .* BF(:,month)  .* ccFineLitter(:,month);
combusted_surfmic       = hsurfmicpool     .* BF(:,month)  .* ccFineLitter(:,month);

% FIRE: the non-combusted parts
nonCombusted_leaf       = hleafpool        .* BF(:,month)  .* (1 - ccLeaf(:,month)) .* mortality_tree;
nonCombusted_froot      = hfrootpool       .* BF(:,month)   * mortality_hfroot;

% FIRE flux from not combusted parts to other pools
hsurfstrpool           = hsurfstrpool + nonCombusted_leaf  .* (1 - metabfract);
hsurfmetpool           = hsurfmetpool + nonCombusted_leaf  .* metabfract;   
hsoilstrpool           = hsoilstrpool + nonCombusted_froot .* (1 - metabfract);
hsoilmetpool           = hsoilmetpool + nonCombusted_froot .* metabfract;   

% FIRE 
hleafpool              = hleafpool        - combusted_leaf        - nonCombusted_leaf;
hfrootpool             = hfrootpool                               - nonCombusted_froot;
hsurfstrpool           = hsurfstrpool     - combusted_surfstr;
hsurfmetpool           = hsurfmetpool     - combusted_surfmet;
hsurfmicpool           = hsurfmicpool     - combusted_surfmic;

% calculate fluxes
if numberOfAgeClasses==1
	hresp          = resppool_surfstr + resppool_surfmet + resppool_surfmic + resppool_soilstr + ...
                                resppool_soilmet + resppool_soilmic + resppool_slow + resppool_armored;
	hcomb          = combusted_leaf + combusted_surfstr + combusted_surfmet + combusted_surfmic;
	hherb          = herbivory;
else
    if ageClass == 1
        hresp = zeros(numberOfGridCells,1); 
        hcomb = zeros(numberOfGridCells,1); 
        hherb = zeros(numberOfGridCells,1); 
        hbiof = zeros(numberOfGridCells,1);
    end
	hresp          = hresp + (resppool_surfstr + resppool_surfmet + resppool_surfmic + resppool_soilstr + resppool_soilmet + ...
                                resppool_soilmic + resppool_slow + resppool_armored) ./ numberOfAgeClasses;
	hcomb          = hcomb + (combusted_leaf + combusted_surfstr + combusted_surfmet + ...
                                combusted_surfmic) ./ numberOfAgeClasses;
	hherb          = hherb + herbivory ./ numberOfAgeClasses;
end    
% *************************************************************************
