% get soil texture parameters
io = soiltext==1; clay(io) = 0.20; silt(io) = 0.20; sand(io) = 0.60;
io = soiltext==2; clay(io) = 0.09; silt(io) = 0.08; sand(io) = 0.83;
io = soiltext==3; clay(io) = 0.20; silt(io) = 0.20; sand(io) = 0.60;
io = soiltext==4; clay(io) = 0.30; silt(io) = 0.33; sand(io) = 0.37;
io = soiltext==5; clay(io) = 0.48; silt(io) = 0.25; sand(io) = 0.27;
io = soiltext==6; clay(io) = 0.67; silt(io) = 0.17; sand(io) = 0.17;
io = soiltext==7; clay(io) = 0.20; silt(io) = 0.20; sand(io) = 0.60;

% get litter characteristic parameters and
% turnover times (leaf and rootage, woodage)
io = veg==1;  litcn(io) = 40; lignin(io) = 0.20; lrage(io) = 1.00; woodage(io) = 75;
io = veg==2;  litcn(io) = 50; lignin(io) = 0.20; lrage(io) = 1.00; woodage(io) = 75;
io = veg==3;  litcn(io) = 65; lignin(io) = 0.22; lrage(io) = 1.00; woodage(io) = 75;
io = veg==4;  litcn(io) = 80; lignin(io) = 0.25; lrage(io) = 3.80; woodage(io) = 75;
io = veg==5;  litcn(io) = 50; lignin(io) = 0.20; lrage(io) = 1.00; woodage(io) = 75;
io = veg==6;  litcn(io) = 50; lignin(io) = 0.15; lrage(io) = 1.00; woodage(io) = 50;
io = veg==7;  litcn(io) = 50; lignin(io) = 0.10; lrage(io) = 1.00; woodage(io) = 50;
io = veg==9;  litcn(io) = 65; lignin(io) = 0.20; lrage(io) = 1.00; woodage(io) = 50;
io = veg==10; litcn(io) = 50; lignin(io) = 0.15; lrage(io) = 1.00; woodage(io) = 50;
io = veg==11; litcn(io) = 50; lignin(io) = 0.15; lrage(io) = 1.00; woodage(io) = 50;
io = veg==12; litcn(io) = 40; lignin(io) = 0.10; lrage(io) = 1.00; woodage(io) = 50;
% Some of these numbers have been CHANGED from original CASA, according to an
% extensive study done for BIOME-BGC (see White, M. A., P. E. Thornton, S.
% W. Running, and R. R. Nemani. 2000. Parameterization and Sensitivity 
% Analysis of the BIOME-BGC Terrestrial Ecosystem Model: Net Primary Production Controls. 
% Earth Interactions 4(3):1-85.)

% here are the original values:
% io = veg==1;  litcn(io) = 40; lignin(io) = 0.20; lrage(io) = 1.80; woodage(io) = 41;
% io = veg==2;  litcn(io) = 50; lignin(io) = 0.20; lrage(io) = 1.20; woodage(io) = 58;
% io = veg==3;  litcn(io) = 65; lignin(io) = 0.22; lrage(io) = 1.20; woodage(io) = 58;
% io = veg==4;  litcn(io) = 80; lignin(io) = 0.25; lrage(io) = 5.00; woodage(io) = 42;
% io = veg==5;  litcn(io) = 50; lignin(io) = 0.20; lrage(io) = 1.80; woodage(io) = 27;
% io = veg==6;  litcn(io) = 50; lignin(io) = 0.15; lrage(io) = 1.80; woodage(io) = 25;
% io = veg==7;  litcn(io) = 50; lignin(io) = 0.10; lrage(io) = 1.50; woodage(io) = 25;
% io = veg==9;  litcn(io) = 65; lignin(io) = 0.20; lrage(io) = 1.00; woodage(io) = 5.5;
% io = veg==10; litcn(io) = 50; lignin(io) = 0.15; lrage(io) = 2.80; woodage(io) = 40;
% io = veg==11; litcn(io) = 50; lignin(io) = 0.15; lrage(io) = 1.00; woodage(io) = 25;
% io = veg==12; litcn(io) = 40; lignin(io) = 0.10; lrage(io) = 1.00; woodage(io) = 25;


% calculate lignin tot nitrogen ratio
ln = (litcn .* lignin) ./ 0.45;

% calculate the rate constants for each pool
annK_leaf       = 1./lrage;
annK_wood       = 1./woodage;
annK_froot      = 1./lrage;
annK_hleaf      = 2;        % turnover time for grass leaves and roots
annK_hfroot     = 2;        % turnover time for grass leaves and roots
annK_surfmet    = 14.8;
annK_surfstr    = 3.9;
annK_soilmet    = 18.5;
annK_soilstr    = 4.8;
annK_cwd        = 0.25;
annK_surfmic    = 6;
annK_soilmic    = 7.3;
annK_slow       = 0.2;
annK_armored    = 0.0045;

% scale the rate constants to their monthly values (not just divided bij 12)
K_wood  = 1-((exp(-annK_wood)).^(1/12));
K_froot = 1-((exp(-annK_froot)).^(1/12));  
K_leaf  = 1-((exp(-annK_leaf)).^(1/12));  

K_hleaf   = 1-((exp(-annK_hleaf))^(1/12));
K_hfroot   = 1-((exp(-annK_hfroot))^(1/12));

K_surfmet = 1-((exp(-annK_surfmet))^(1/12));
K_surfstr = 1-((exp(-annK_surfstr))^(1/12));
K_soilmet = 1-((exp(-annK_soilmet))^(1/12));
K_soilstr = 1-((exp(-annK_soilstr))^(1/12));
K_cwd     = 1-((exp(-annK_cwd))^(1/12));    
K_surfmic = 1-((exp(-annK_surfmic))^(1/12));
K_soilmic = 1-((exp(-annK_soilmic))^(1/12));
K_slow    = 1-((exp(-annK_slow))^(1/12));
K_armored = 1-((exp(-annK_armored))^(1/12));

% microbial afficiencies for particular flows
eff_surfstr2slow    = 0.7;
eff_surfstr2surfmic = 0.4;
eff_soilstr2slow    = 0.7;
eff_soilstr2soilmic = 0.45;
eff_cwd2slow        = 0.7;
eff_cwd2surfmic     = 0.4;
eff_surfmic2slow    = 0.4;
eff_surfmet2surfmic = 0.4;
eff_soilmet2soilmic = 0.45;
eff_slow2soilmic    = 0.45;
eff_armored2soilmic = 0.45;
eff_soilmic2slow    = 0.85 - (0.68.*(silt + clay));

% determine what fraction of the litter will be metabolic
metabfract = 0.85 - (0.018 .* ln);
metabfract(metabfract<0) = 0;

% get the fraction of the carbon in the structural litter pools
% that will be from lignin
structuralLignin = ((lignin .* 0.65) ./ 0.45) ./ (1 - metabfract);

% estimate of the lignin content of wood carbon is 40%
woodligninfract  = 0.4;
lignineffect     = exp(-3 .* structuralLignin);

% calculate vegetation-dependent factors used in the belowground module
soilmicDecayFactor = (1 - (0.75 .* (silt + clay)));
slowDecayFactor    = ones(numberOfGridCells,1);  
armoredDecayFactor = ones(numberOfGridCells,1);

% decay is faster in agriculture gridcells (till etc)
io = veg==12;
soilmicDecayFactor(io) = soilmicDecayFactor(io) .* 1.25;
slowDecayFactor(io)    = slowDecayFactor(io)    .* 1.5;
armoredDecayFactor(io) = armoredDecayFactor(io) .* 1.5;

fid = 0.003 - (0.009 .* clay);
decayClayFactor = fid;
decayClayFactor(fid < 0) = 0;

