diff --git a/README.md b/README.md index 45e330d..b3fdd8d 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,25 @@ # ed-allometry-unittest This is a matlab unit-testing framework for unit-testing various allometry schemes in the ED model. These are offline (not coupled to the CLM-ED, ED2 or other) tests. + +Contents: + +allom_lib_v3.m: The library that is mirrored in moderen fortran inside the model. + +drive_allomtests.m: The user invokes tests through this script. This is the driver script and contains some control parameters. + +gen_param_instance.m: This is the control card that defines the parameters associated with the different schemes. + +The following are various plotting functions that could un-doubtedly be re-written more efficiently into one script with variable controls: + +functions/plot_multicase_bag.m +functions/plot_multicase_blmax.m +functions/plot_multicase_blmax_o_dbagdh.m +functions/plot_multicase_blmaxdi.m +functions/plot_multicase_bsap.m +functions/plot_multicase_bsapbag.m +functions/plot_multicase_bsapid.m +functions/plot_multicase_dbhe.m +functions/plot_multicase_h.m +functions/plot_singlecase_cfractions.m + +gen_param_instance.m diff --git a/allom_lib_v3.m b/allom_lib_v3.m new file mode 100644 index 0000000..75afeb5 --- /dev/null +++ b/allom_lib_v3.m @@ -0,0 +1,1056 @@ + +%========================================================================== +% +% allom_lib_v3.m +% A library of functions that calculate plant allometry and their +% derivatives. All relationships are related to diameter, all +% derivatives are assumed to be with respect to change in diameter with +% units of [cm]. In some cases multiple areguments will be provided, yet +% those arguments in those cases are trivially determined from diameter. +% +% Each function presented in this library is written to also return the +% derivative with respect to diameter using a logical switch "dswitch" +% (derivative-switch). With one exception, the h2d function returns the +% change in diameter with respect to height. +% +% The name convention of the functions follows the form d2... Which +% indicates "diameter to ...". Allometries for the following variables are +% calculated: +% h: height [m] +% bag: biomass above ground [kgC] (aka AGB) +% blmax: biomass in leaves when leaves are "on allometry" +% this also is the "pre-trimmed" value, which is the maximum +% or potential leaf mass a tree may have [kgC] +% bcr: biomass in coarse roots [kgC] (belowground sap+dead wood, no fines) +% bfrmax: biomass in fine roots when "on allometry" [kgC] +% bsap: biomass in sapwood (above and below) [kgC] +% bdead: biomass (above and below) in the structural pool [kgC] +% +% "on allometry" assumes the plant has leaves flushed, and has had +% sufficient carbon to meet maintenance turnover. +% +% The following model traits are used: +% h_mode, string, height allometry function type +% blmax_mode, string, maximum leaf allometry function type +% bfrmax_mode, string, maximum root allometry function type +% bag_mode, string, AGB allometry function type +% bcr_mode, string, coarse root allometry function type +% bsap_mode, string, sapwood allometry function type +% wood_density, real, mean stem wood specific gravity (heart,sap,bark) +% latosa_int, real, leaf area to sap area ratio, intercept [m2/cm2] +% latosa_slp, real, leaf area to sap area ratio, slope on diameter +% [m2/cm2/cm] +% c2b, real, carbon to biomass ratio (~2.0) +% d2h1, real, parameter 1 for d2h allometry (intercept) +% d2h2, real, parameter 2 for d2h allometry (slope) +% d2h3, real, parameter 3 for d2h allometry (optional) +% eclim, real, climatological influence parameter for d2h allometry +% bl_min, real, leaf biomass of smallest tracked plant [kgC] +% froot_leaf, real, fine root biomass per leaf biomass ratio [kgC/kgC] +% ag_biomass, real, the fraction of stem above ground [-] +% dbh_adult, real, the lower dbh bound where adult allometry starts (10m) +% h_max, real, maximum height of a functional type/group +% d2bl1, real, parameter 1 for d2bl allometry (intercept) +% d2bl2, real, parameter 2 for d2bl allometry (slope) +% d2bl3, real, parameter 3 for d2bl allometry (optional) +% h_min, real, the height associated with newly recruited plant [m] +% dbh_min, real, the dbh associated with a newly recruited plant [cm] +% dbh_max, real, the diameter associated with maximum height [cm] +% diagnosed from hmax using non-asymptotic functions +% +%========================================================================== + +function [p_h,p_bag,p_blmax,p_h2d, ... + p_bsap,p_bcr,p_bfrmax,p_bdead] = allom_lib_v3 + +p_h = @f_h; +p_bag = @f_bag; +p_blmax = @f_blmax; +p_h2d = @f_h2d; +p_bcr = @f_bcr; +p_bsap = @f_bsap; +p_bfrmax = @f_bfrmax; +p_bdead = @f_bdead; +end + +% ========================================================================= +% Parameter Checks and Defaults (subroutine) +% ========================================================================= + + + + + +% ========================================================================= +% Generic height to diameter interface +% ========================================================================= + +function [d,dddh] = f_h2d(h,traitp,ipf) +switch traitp.h_mode{ipf} + case {'chave14'} + [d,dddh] = h2d_chave2014(h,traitp,ipf); + case {'poorter06'} + [d,dddh] = h2d_poorter2006(h,traitp,ipf); + case {'2par_pwr'} + [d,dddh] = h2d_2pwr(h,traitp,ipf); + case {'obrien'} + [d,dddh] = h2d_obrien(h,traitp,ipf); + otherwise + display('Unknown D-2-H Allometry'); + pause; + return; +end +end + +% ========================================================================= +% Generic height interface +% ========================================================================= + +function [h,dhdd] = f_h(d,traitp,ipf) +switch traitp.h_mode{ipf} + case {'chave14'} + [h,dhdd] = d2h_chave2014(d,traitp,ipf); + case {'poorter06'} + [h,dhdd] = d2h_poorter2006(d,traitp,ipf); + case {'2par_pwr'} + [h,dhdd] = d2h_2pwr(d,traitp,ipf); + case {'obrien'} + [h,dhdd] = d2h_obrien(d,traitp,ipf); + otherwise + display('Unknown H Allometry'); + pause; + return; +end +end + +% ========================================================================= +% Generic AGB interface +% ========================================================================= + +function [bag,dbagdd] = f_bag(d,h,traitp,ipf) +switch traitp.bag_mode{ipf} + case {'chave14'} + [bag,dbagdd] = dh2bag_chave2014(d,h,traitp,ipf); + case {'2par_pwr'} + % Switch for woodland dbh->drc + [bag,dbagdd] = d2bag_2pwr(d,traitp,ipf); + case {'salda'} + [bag,dbagdd] = dh2bag_salda(d,h,traitp,ipf); + otherwise + display('Unknown BAG Allometry'); + pause; + return; +end +end + +% ========================================================================= +% Generic diameter to maximum leaf biomass interface +% ========================================================================= + +function [blmax,dblmaxdd] = f_blmax(d,h,traitp,ipf) + +switch traitp.blmax_mode{ipf} + case{'salda'} + [blmax,dblmaxdd] = d2blmax_salda(d,traitp,ipf); + case{'2par_pwr'} + [blmax,dblmaxdd] = d2blmax_2pwr(d,traitp,ipf); + case{'2par_pwr_asym'} + [blmax,dblmaxdd] = dh2blmax_2pwr_spline_asym(d,h,traitp,ipf); + case{'2par_pwr_hcap'} + [blmax,dblmaxdd] = d2blmax_2pwr_hcap(d,traitp,ipf); + otherwise + display(sprintf('Uknown blmax Allometry: %s',traitp.blmax_mod{ipf})); + pause; + return; +end +end + +% ========================================================================= +% Generic sapwood biomass interface +% ========================================================================= + +function [bsap,dbsapdd] = f_bsap(d,h,blmax,dblmaxdd,dhdd,traitp,ipf) + +switch traitp.bsap_mode{ipf} + % --------------------------------------------------------------------- + % Currently both sapwood area proportionality methods use the same + % machinery. The only differences are related to the parameter + % checking at the beginning. For constant proportionality, the slope + % of the la:sa to diameter line is zero. + % --------------------------------------------------------------------- + case{'constant','dlinear'} + [bsap,dbsapdd] = bsap_dlinear(d,h,blmax,dblmaxdd,dhdd,traitp,ipf); + otherwise + display('Uknown BSAP Allometry'); + pause; + return; +end +end + +% ========================================================================= +% Generic coarse root biomass interface +% ========================================================================= + +function [bcr,dbcrdd] = f_bcr(d,bag,dbagdd,traitp,ipf) + +switch traitp.bcr_mode{ipf} + case{'constant'} + [bcr,dbcrdd] = bcr_const(d,bag,dbagdd,traitp,ipf); + otherwise + display('Uknown D-2-BCR Allometry'); + pause; + return; +end +end + +% ========================================================================= +% Generic maximum fine root biomass interface +% ========================================================================= + +function [bfrmax,dbfrmaxdd] = f_bfrmax(d,blmax,dblmaxdd,traitp,ipf) + +switch traitp.bfrmax_mode{ipf} + case{'constant'} + [bfrmax,dbfrmaxdd] = bfrmax_const(d,blmax,dblmaxdd,traitp,ipf); + otherwise + display('Uknown D-2-BFR Allometry'); + pause; + return; +end +end + +% ========================================================================= +% Dead biomass interface +% ========================================================================= + +function [bdead,dbdeaddd] = f_bdead(bag,bcr,blmax,bsap, ... + dbagdd,dbcrdd,dblmaxdd,dbsapdd) %#ok + +% bdead is diagnosed as the mass balance from all other pools +% and therefore, no options are necessary +% We are ignoring blmax right now, because it is insignificant in large +% trees and may cause negatives in treelets and saplings + +%bdead = bag+bcr-blmax-bsap; +bdead = bag+bcr-bsap; +dbdeaddd = dbagdd+dbcrdd-dbsapdd; + +end + +% ========================================================================= +% Specific bfrmax relationships +% ========================================================================= + +function [bfrmax,dbfrmaxdd] = bfrmax_const(d,blmax,dblmaxdd,traitp,ipf) %#ok + +froot_leaf = traitp.froot_leaf(ipf); % fine root biomass [kgC] per leaf biomass [kgC] + +bfrmax = blmax.*froot_leaf; +% dbfr/dd = dbfrmax/dblmax * dblmax/dd +dbfrmaxdd = froot_leaf.*dblmaxdd; + +end + + +% ========================================================================= +% Specific bcr relationships +% ========================================================================= + +function [bcr,dbcrdd] = bcr_const(d,bag,dbagdd,traitp,ipf) %#ok + +ag_biomass = traitp.ag_biomass(ipf); % Above-ground biomass fraction + +% btot = bag + bcr; +% bag = btot*ag_biomass; +% bag/ag_biomass = bag + bcr; +% bcr = bag*(1/ag_biomass-1) + +bcr = bag*(1/ag_biomass-1); +% Derivative +% dbcr/dd = dbcr/dbag * dbag/dd +dbcrdd = (1/ag_biomass-1).*dbagdd; + +end + + +% ========================================================================= +% Specific d2bsap relationships +% ========================================================================= + + +function [bsap,dbsapdd] = bsap_dlinear(dbh,h,blmax,dblmaxdd,dhdd,traitp,ipf) + +% ------------------------------------------------------------------------- +% Calculate sapwood biomass based on leaf area to sapwood area +% proportionality. In this function, the leaftosapwood area is a function +% of plant size, see Calvo-Alvarado and Bradley Christoferson +% In this case: parameter latosa (from constant proportionality) +% is the intercept of the diameter function. +% +% For very small plants, the fraction can get very large, so cap the amount +% of sapwood at X% of agb-bleaf; +% ------------------------------------------------------------------------- + +max_agbfrac = 0.75; + +gtokg = 1000.0; % gram to kg conversion +cm2tom2 = 10000.0; % cm^2 to m^2 conversion +mg2kg = 1000.0; % megagrams to kilograms conversion [kg/Mg] +latosa_int = traitp.latosa_int(ipf); % leaf area to sap area ratio [m2/cm2] +latosa_slp = traitp.latosa_slp(ipf); % slope of diameter to latosa line [m2/cm2/cm] +sla = traitp.slatop(ipf); % Specific leaf area [m2/gC] +wood_density = traitp.wood_density(ipf); % Wood density (specific grav) [Mg/m3] +c2b = traitp.c2b(ipf); % carbon to biomass ratio + +% ------------------------------------------------------------------------- +% Calculate sapwood biomass per linear height and kgC of leaf [m-1] +% Units: (1/latosa)* slatop* gtokg * cm2tom2 / c2b * mg2kg * dens +% [cm2/m2] *[m2/gC]*[1000gC/1kgC]*[1m2/10000cm2] /[kg/kgC]*[kg/Mg]*[Mg/m3] +% ->[cm2/gC] +% ->[cm2/kgC] +% ->[m2/kgC] +% ->[m2/kg] +% ->[m2/Mg] +% ->[/m] +% ------------------------------------------------------------------------- + +latosa = latosa_int + dbh*latosa_slp; + +hbl2bsap = sla*gtokg*wood_density*mg2kg/(latosa*c2b*cm2tom2); +[bag,dbagdd] = f_bag(dbh,h,traitp,ipf); + +bsap = min(max_agbfrac*bag,hbl2bsap * h * blmax); + +% Derivative +% dbldmaxdd is deriv of blmax wrt dbh (use directives to check oop) +% dhdd is deriv of height wrt dbh (use directives to check oop) + +if(bsap=dbh_sap && dbh=dbh_adult && dbh=dbh_adult) + blmax = d2bl1.*dbh.^d2bl2 ./ c2b; + dblmaxdd = d2bl1.*d2bl2.*dbh.^(d2bl2-1.0) ./ c2b; +else + bleaf_ad = d2bl1.*dbh_adult.^d2bl2 ./c2b; + a2_small = log(bleaf_ad./bl_min)./log(dbh_adult./dbh_min); + a1_small = bleaf_ad.*c2b/(dbh_adult.^a2_small); + blmax = a1_small.*dbh.^a2_small ./ c2b; + dblmaxdd = a1_small.*a2_small.*dbh.^(a2_small-1.0) ./ c2b; +end +end + +function [blmax,dblmaxdd] = dh2blmax_2pwr_spline_asym(dbh,h,traitp,ipf) + +d2bl1 = traitp.d2bl1(ipf); +d2bl2 = traitp.d2bl2(ipf); +c2b = traitp.c2b(ipf); +dbh_adult = traitp.dbh_adult(ipf); + +dbh_sap = traitp.dbh_sap(ipf); +d2bl1_sap = traitp.d2bl1_sap(ipf); %0.0201; +d2bl2_sap = traitp.d2bl2_sap(ipf); %3.1791; + +[dbh_eff,ddbhedh] = f_h2d(h,traitp,ipf); + +% In this version of the 2 parameter power fit of diameter to maximum leaf +% biomass; leaf biomass is capped based on maximum height. However, leaf +% allometry is based on diameter. So the actual asymptoted height +% calculates the diameter that would had generated that height in a +% condition where no asymptote had occured. This is called dbh_eff. + +% There are two portions of the curve. Trees that are adult stature and +% larger and thus applicable to the published regressions, and those +% smaller than that are thus applicable to an extrapolation to known leaf +% biomass found in saplings. + + +if(dbh_eff=dbh_sap && dbh_eff=dbh_adult && dbh<=dbh_maxh) + blmax = d2bl1.*dbh.^d2bl2./c2b; + dblmaxdd = d2bl2.*d2bl1.*dbh.^(d2bl2-1)./c2b; +elseif(dbh>dbh_maxh) + blmax = d2bl1.*dbh_maxh.^d2bl2./c2b; + dblmaxdd = 0.0; +else + bleaf_ad = d2bl1.*dbh_adult.^d2bl2 ./c2b; + d2bl2_small = log(bleaf_ad./bl_min)./log(dbh_adult./dbh_min); + d2bl1_small = bleaf_ad.*c2b/(dbh_adult.^d2bl2_small); + blmax = d2bl1_small.*dbh.^d2bl2_small ./ c2b; + dblmaxdd = d2bl1_small.*d2bl2_small.*dbh.^(d2bl2_small-1)./ c2b; +end +end + +% ========================================================================= +% Diameter to height (D2H) functions +% ========================================================================= + +function [h,dhdd] = d2h_chave2014(dbh,traitp,ipf) + +% "d2h_chave2014" +% "dbh to height via Chave et al. 2014" + +% This function calculates tree height based on tree diameter and the +% environmental stress factor "E", as per Chave et al. 2015 GCB +% As opposed to previous allometric models in ED, in this formulation +% we do not impose a hard cap on tree height. But, maximum_height +% is an important parameter, but instead of imposing a hard limit, in +% the new methodology, it will be used to trigger a change in carbon +% balance accounting. Such that a tree that hits its maximum height will +% begin to route available NPP into seed and defense respiration. +% +% The stress function is based on the geographic location of the site. If +% a user decides to use Chave2015 allometry, the E factor will be read in +% from a global gridded dataset and assigned for each ED patch (note it +% will be the same for each ED patch, but this distinction will help in +% porting ED into different models (patches are pure ED). It +% assumes that the site is within the pan-tropics, and is a linear function +% of climatic water deficit, temperature seasonality and precipitation +% seasonality. See equation 6b of Chave et al. +% The relevant equation for height in this function is 6a of the same +% manuscript, and is intended to pair with diameter to relate with +% structural biomass as per equation 7 (in which H is implicit). +% +% Chave et al. Improved allometric models to estimate the abovegroud +% biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2015. +% +% ========================================================================= + +%eclim: Chave's climatological influence parameter "E" + +d2h1 = traitp.d2h1(ipf); % (alias) +d2h2 = traitp.d2h2(ipf); % (alias) +d2h3 = traitp.d2h3(ipf); % (alias) +eclim = traitp.eclim(ipf); % (alias) +dbh_maxh = traitp.dbh_maxh(ipf); +ddbh = 0.1; % 1-mm + +% For the non derivative solution, if the tree is large and +% close to any cap that is imposed, then we need to perform a +% step-integration because the asymptotic function makes the indefinite +% integral incredibly messy. Thus we use an Euler step, yes ugly, +% but it is a simple function so don't over-think it (RGK). + +k = 0.25; + +if(dbh>0.5*dbh_maxh) + + dbh0=0.5*dbh_maxh; + h = exp( d2h1 - eclim + d2h2*log(dbh0) + d2h3*log(dbh0).^2.0 ); + while(dbh00.5*dbh_maxh) + dbh0=0.5*dbh_maxh; + h = d2h1.*dbh0.^d2h2; + while(dbh0=dbh_sap && dbh=dbh_hmax) + h = 10.^(log10(dbh_hmax)*d2h1+d2h2); + dhdd = 0; +else + h = 10.^(log10(dbh)*d2h1+d2h2); + dhdd = d2h1.*10.^d2h2.*dbh.^(d2h1-1.0); +end + +end + + + +% ========================================================================= +% Diameter 2 above-ground biomass +% ========================================================================= + +function [bag,dbagdd] = dh2bag_chave2014(dbh,h,traitp,ipf) + +% ========================================================================= +% This function calculates tree structural biomass from tree diameter, +% height and wood density. +% +% Chave et al. Improved allometric models to estimate the abovegroud +% biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2015. +% +% Input arguments: +% dbh: Diameter at breast height [cm] +% rho: wood specific gravity (dry wood mass per green volume) +% height: total tree height [m] +% a1: structural biomass allometry parameter 1 (intercept) +% a2: structural biomass allometry parameter 2 (slope) +% Output: +% bag: Total above ground biomass [kgC] +% +% ========================================================================= + +d2bag1 = traitp.d2bag1(ipf); +d2bag2 = traitp.d2bag2(ipf); +wood_density = traitp.wood_density(ipf); +c2b = traitp.c2b(ipf); + +bag = (d2bag1 * (wood_density*dbh.^2*h).^d2bag2)/c2b; + +[~,dhdd] = f_h(dbh,traitp,ipf); +dbagdd1 = (d2bag1.*wood_density.^d2bag2)./c2b; +dbagdd2 = d2bag2.*dbh.^(2.*d2bag2).*h.^(d2bag2-1.0).*dhdd; +dbagdd3 = h.^d2bag2.*2.*d2bag2.*dbh.^(2.*d2bag2-1); +dbagdd = dbagdd1.*(dbagdd2 + dbagdd3); + +end + + +function [bag,dbagdd] = d2bag_2pwr(diam,traitp,ipf) + +% ========================================================================= +% This function calculates tree above ground biomass according to 2 +% parameter power functions. (slope and intercepts of a log-log +% diameter-agb fit: +% +% These relationships are typical for temperate/boreal plants in North +% America. Parameters are available from Chojnacky 2014 and Jenkins 2003 +% +% Note that we are using an effective diameter here, as Chojnacky 2014 +% and Jenkins use diameter at the base of the plant for "woodland" species +% The diameters should be converted prior to this routine if drc. +% +% Input arguments: +% diam: effective diameter (dbh or drc) in cm +% FOR SPECIES THAT EXPECT DCM, THIS NEEDS TO BE PRE-CALCULATED!!!! +% Output: +% agb: Total above ground biomass [kgC] +% +% ========================================================================= +% Aceraceae, Betulaceae, Fagaceae and Salicaceae comprised nearly +% three-fourths of the hardwood species (Table 3) +% +% Fabaceae and Juglandaceae had specific gravities .0.60 and were +% combined, as were Hippocastanaceae and Tilaceae with specific gravities +% near 0.30. The remaining 9 families, which included mostly species with +% specific gravity 0.45–0.55, were initially grouped to construct a general +% hardwood taxon for those families having few published biomass equa- +% tions; however, 3 warranted separation, leaving 6 families for the general +% taxon. +% bag = exp(b0 + b1*ln(diameter))/c2b; +% ========================================================================= + +d2bag1 = traitp.d2bag1(ipf); +d2bag2 = traitp.d2bag2(ipf); +c2b = traitp.c2b(ipf); + +%max_dbh = traitp.maxdbh(ipf); +%if(diam>1.10*max_dbh) +% display('-----------------------------------------------------'); +% display('Tree diameter is 10% larger than diameter where height'); +% display('hits maximum. However, you specified an AGB allometry'); +% display('that does not assume capping. Please consider '); +% display('re-evaluating your allometric assumptions, growth'); +% display('formulations or maximum height'); +% display('------------------------------------------------------'); +%end + +bag = (d2bag1 * diam.^d2bag2)./c2b; +dbagdd = (d2bag2.*d2bag1.*diam.^(d2bag2-1.0))./c2b; + +end + +function [bag,dbagdd] = dh2bag_salda(dbh,h,traitp,ipf) + +% In the interest of reducing the number of model parameters, and since +% Saldarriaga 1988 seems as though it is being deprecated, we will use +% hard-wired parameter for dh2bag_salda, which would had required 4 +% variable parameters. + +d2bag1 = 0.06896; +d2bag2 = 0.572; +d2bag3 = 1.94; +d2bag4 = 0.931; +c2b = traitp.c2b(ipf); +wood_density = traitp.wood_density(ipf); + +bag = d2bag1*(h^d2bag2)*(dbh^d2bag3)*(wood_density.^d2bag4)./c2b; + +% bag = a1 * h^a2 * d^a3 * r^a4 +% dbag/dd = a1*r^a4 * d/dd (h^a2*d^a3) +% dbag/dd = a1*r^a4 * [ d^a3 *d/dd(h^a2) + h^a2*d/dd(d^a3) ] +% dbag/dd = a1*r^a4 * [ d^a3 * a2*h^(a2-1)dh/dd + h^a2*a3*d^(a3-1)] + +term1 = d2bag1*(wood_density.^d2bag4)./c2b; +term2 = (h^d2bag2)*d2bag3*dbh^(d2bag3-1.0); +[~,dhdd] = f_h(dbh,traitp,ipf); +term3 = d2bag2*h^(d2bag2-1)*(dbh^d2bag3)*dhdd; +dbagdd = term1*(term2+term3); + +end + +function [dbhe,ddbhedh] = h2d_chave2014(h,traitp,ipf) + +d2h1 = traitp.d2h1(ipf); % (alias) +d2h2 = traitp.d2h2(ipf); % (alias) +d2h3 = traitp.d2h3(ipf); % (alias) +eclim = traitp.eclim(ipf); % (alias) +ar = d2h1-eclim; + +eroot = (-d2h2 + sqrt(d2h2^2 + 4*log(h.*exp(-ar))*d2h3))./(2*d2h3); +dbhe = exp(eroot); + +% How about just inverting the derivative at phi? without asyms +eroot = (-d2h2 + sqrt(d2h2^2 + 4*log(h.*exp(-ar))*d2h3))./(2*d2h3); +dbh1 = exp(eroot); +dhpdd = exp(ar).*( d2h3.*2.0.*dbh1.^(d2h2-1.0).*log(dbh1).* ... + exp(d2h3.*log(dbh1).^2) + d2h2.*dbh1.^(d2h2-1.0).* ... + exp(d2h3.*log(dbh1).^2.0) ); + +ddbhedh = 1/dhpdd; + +% term1 = exp(-d2h2./(2*d2h3)); +% term2 = exp(d2h2.^2./(4.*d2h3.^2)); +% term3 = exp(-ar/d2h3); +% term4 = h.^(1/d2h3-1.0)./(d2h3); +% dbh = term1*term2*term3*term4; +end + + +function [dbh,ddbhdh] = h2d_poorter2006(h,traitp,ipf) + +% ------------------------------------------------------------------------- +% Note that the height to diameter conversion in poorter is only necessary +% when initializing saplings. In other methods, height to diameter is +% useful when defining the dbh at which point to asymptote, as maximum +% height is the user set parameter. This function should not need to set a +% dbh_max parameter for instance, but it may end up doing so anyway, even +% if it is not used, do to poor filtering. The poorter et al. d2h and h2d +% functions are already asymptotic, and the parameter governing maximum +% height is the d2h1 parameter. Note as dbh gets very large, the +% exponential goes to zero, and the maximum height approaches d2h1. +% However, the asymptote has a much different shape than the logistic, so +% therefore in the Poorter et al functions, we do not set d2h1 == h_max. +% That being said, if an h_max that is greater than d2h1 is passed to this +% function, it will return a complex number. During parameter +% initialization, a check will be placed that forces: +% h_max = d2h1*0.98 +% ------------------------------------------------------------------------- + + +d2h1 = traitp.d2h1(ipf); %(alias) +d2h2 = traitp.d2h2(ipf); %(alias) +d2h3 = traitp.d2h3(ipf); %(alias) + +% ------------------------------------------------------------------------- +% h = a1.*(1 - exp(a2.*dbh.^a3)) +% h = a1 - a1*exp(a2*dbh^a3) +% a1-h = a1*exp(a2*dbh^a3) +% (a1-h)/a1 = exp(a2*dbh^a3) +% log(1-h/a1) = a2*dbh^a3 +% [log(1-h/a1)/a2]^(1/a3) = dbh +% +% derivative dd/dh +% dd/dh = [log((a1-h)/a1)/a2]^(1/a3)' +% = (1/a3)*[log((a1-h)/a1)/a2]^(1/a3-1)* [(log(a1-h)-log(a1))/a2]' +% = (1/a3)*[log((a1-h)/a1)/a2]^(1/a3-1) * (1/(a2*(h-a1)) +% dd/dh = -((log(1-h/a1)/a2).^(1/a3-1))/(a2*a3*(a1-h)) +% ------------------------------------------------------------------------- + +dbh = (log(1.0-h/d2h1)./d2h2).^(1.0/d2h3); +ddbhdh = -((log(1-h/d2h1)/d2h2).^(1/d2h3-1))./(d2h2*d2h3*(d2h1-h)); +end + + +function [dbh,ddbhdh] = h2d_2pwr(h,traitp,ipf) + +d2h1 = traitp.d2h1(ipf); +d2h2 = traitp.d2h2(ipf); + +%h = a1.*dbh.^a2; +dbh = (h/d2h1).^(1/d2h2); +% dbh = (1/a1).^(1/a2).*h^(1/a2); +ddbhdh = (1/d2h2).*(1/d2h1).^(1/d2h2).*h^(1/d2h2-1.0); +end + + +function [dbh,ddbhdh] = h2d_obrien(h,traitp,ipf) + +d2h1 = traitp.d2h1(ipf); +d2h2 = traitp.d2h2(ipf); +h_max = traitp.h_max(ipf); +dbh_sap = traitp.dbh_sap(ipf); +dbh_adult = traitp.dbh_adult(ipf); +d2h1_sap = traitp.d2h1_sap(ipf); %1.0729; +d2h2_sap = traitp.d2h2_sap(ipf); %1.4925; + +% ========================================================================= +% From King et al. 1990 at BCI for saplings +% log(d) = a + b*log(h) +% d = exp(a) * h^b +% h = (1/exp(a)) * d^(1/b) +% h = d2h1*d^d2h2 where d2h1 = 1/exp(a) d2h2 = 1/b +% d = (h/d2h1)^(1/d2h2); +% For T. tuberculata (canopy tree) a1 = -0.0704, b1 = 0.67 +% ========================================================================= + +h_sap = d2h1_sap * dbh_sap^d2h2_sap; +h_adult = 10.0^(log10(dbh_adult) * d2h1 + d2h2); + +if(h=h_sap && h dbh_hmax +% 5) Sanity Checks +% 6) 0.1% bag < blmax < 10% bag | dbh_min -> dbh_hmax +% 7) 0.1% bag < bfrmax < 10% bag | dbh_min -> dbh_hmax +% 8) 3% bag < bcr < 90% bag | dbh_min -> dbh_hmax +% 9) 0.1% bag < bsap < 60% bag | dbh_min -> dbh_hmax +% 10) Visual tests, plot all functional relationships x-y +% all cases are displayed in each x/y graph for: +% d2h, h2d, d2bag, d2blmax, d2sap, d2bcr, d2bfrmax, d2bdead +% +% ========================================================================= + +% Calculate Results from Cases +% ========================================================================= +for ic=1:nc + + % Initialize Both Integrated and Absolute Quantities + % --------------------------------------------------------------------- + d0 = dbh(ic,1); + [h0,dhdd0] = f_h(d0,traitp,ic); + % Height Integrated + hi(ic,1) = h0; + % AGB + [bagi(ic,1),~] = f_bag(d0,h0,traitp,ic); + % AGB + [bagd(ic,1),dbagd0] = f_bag(d0,h0,traitp,ic); + % blmax (integrated) + [blmaxi(ic,1),~] = f_blmax(d0,h0,traitp,ic); + % blmax (integrated) + [blmaxd(ic,1),~] = f_blmax(d0,0,traitp,ic); + % bfrmax + [bfrmax(ic,1),~] = f_bfrmax(d0,blmaxi(ic,1),0,traitp,ic); + % bcr + [bcr(ic,1),~] = f_bcr(d0,bagi(ic,1),0,traitp,ic); + % bsap (integrated) + [bsapi(ic,1),~] = f_bsap(d0,hi(ic,1),blmaxi(ic,1),0,0,traitp,ic); + % bsap (direct) + [bsapd(ic,1),~] = f_bsap(d0,hi(ic,1),blmaxi(ic,1),0,0,traitp,ic); + % bdead + [bdead(ic,1),~] = f_bdead(bagi(ic,1),bcr(ic,1),blmaxi(ic,1),bsapi(ic,1),0,0,0,0); + % Reverse (effective diameter) + dbhe(ic,1) = dbh(ic,1); + + + blmax_o_dbagdh(ic,1) = blmaxi(ic,1)./(dbagd0/dhdd0); + + % Walk through a range of diameters + for id=2:ndbh + dp = dbh(ic,id-1); + dc = dbh(ic,id); + dd = dc-dp; + + % Height + [~,dhdd] = f_h(dp,traitp,ic); + hi(ic,id) = hi(ic,id-1) + dhdd*dd; + + % Effective diameter( hd2(h) ) + %if(ic==3) + [dbhe(ic,id),~] = f_h2d(hi(ic,id),traitp,ic); + %end + + % AGB (integrated) + [~,dbagdd] = f_bag(dp,hi(ic,id-1),traitp,ic); + bagi(ic,id) = bagi(ic,id-1) + dbagdd*dd; + + % AGB (direct) + [bagd(ic,id),dbagdd] = f_bag(dc,hi(ic,id),traitp,ic); + + % blmax (integrated) + [~,dblmaxdd] = f_blmax(dp,hi(ic,id-1),traitp,ic); + blmaxi(ic,id) = blmaxi(ic,id-1) + dblmaxdd*dd; + + % blmax (direct) + [blmaxd(ic,id),~] = f_blmax(dc,hi(ic,id),traitp,ic); + + % bfrmax + [~,dbfrmaxdd] = f_bfrmax(dp,blmaxi(ic,id-1),dblmaxdd,traitp,ic); + bfrmax(ic,id) = bfrmax(ic,id-1) + dbfrmaxdd*dd; + + % bcr + [~,dbcrdd] = f_bcr(dp,bagi(ic,id-1),dbagdd,traitp,ic); + bcr(ic,id) = bcr(ic,id-1) + dbcrdd*dd; + + % bsap (integrated) + [~,dbsapdd] = f_bsap(dp,hi(ic,id-1),blmaxi(ic,id-1),dblmaxdd,dhdd,traitp,ic); + bsapi(ic,id) = bsapi(ic,id-1) + dbsapdd*dd; + + % bsap (direct) + [bsapd(ic,id),~] = f_bsap(dp,hi(ic,id),blmaxi(ic,id),0,0,traitp,ic); + + if(hi(ic,id)>=traitp.h_max(ic)) + blmax_o_dbagdh(ic,id) = NaN; + else + blmax_o_dbagdh(ic,id) = blmaxi(ic,id)./(dbagdd./dhdd); + end + + %display([bagi(ic,id)+bcr(ic,id)-blmaxi(ic,id)-bsapi(ic,id),bagi(ic,id),bcr(ic,id),blmaxi(ic,id),bsapi(ic,id)]) + %pause; + + % bdead + [~,dbdeaddd] = f_bdead(bagi(ic,id-1),bcr(ic,id-1), ... + blmaxi(ic,id-1),bsapi(ic,id-1), ... + dbagdd,dbcrdd,dblmaxdd,dbsapdd); + bdead(ic,id) = bdead(ic,id-1) + dbdeaddd*dd; + + end + + +end + + + +% Test 1: Make direct plots of each +fig=1; plot_multicase_h(dbh,hi,traitp,fig); +fig=fig+1; plot_multicase_bag(dbh,bagi,traitp,fig); +fig=fig+1; plot_multicase_blmax(dbh,blmaxi,traitp,fig); +%fig=fig+1; plot_multicase_bfrmax(dbh,bfrmaxi,traitp,fig); +%fig=fig+1; plot_multicase_bcr(state,traitp,fig); +fig=fig+1; plot_multicase_bsap(dbh,bsapd,traitp,fig); +fig=fig+1; plot_multicase_bsapid(bsapd,bsapi,traitp,fig); +fig=fig+1; plot_multicase_dbhe(dbh,dbhe,traitp,fig); +fig=fig+1; plot_multicase_blmaxdi(blmaxi,blmaxd,traitp,fig); +fig=fig+1; plot_multicase_bsapbag(dbh,bsapi./bagi,traitp,fig); +fig=fig+1; plot_multicase_blmax_o_dbagdh(dbh,blmax_o_dbagdh,traitp,fig); + +for ic=1:nc + fig=fig+1; plot_singlecase_cfractions(dbh(ic,:),blmaxi(ic,:), ... + bfrmax(ic,:),bcr(ic,:),bsapi(ic,:),bdead(ic,:), ... + traitp,ic,fig) + pause; +end diff --git a/functions/plot_multicase_bag.m b/functions/plot_multicase_bag.m new file mode 100644 index 0000000..0db7aea --- /dev/null +++ b/functions/plot_multicase_bag.m @@ -0,0 +1,20 @@ +function plot_multicase_bag(dbh,bag,traitp,fign) + +[ncases,~] = size(dbh); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(dbh',bag'./1000,'LineWidth',2); +xlabel('DBH [cm]'); +ylabel('Aboveground Biomass [Mg]'); +title(sprintf('Allometry Visual Tests\nIntercomparison of AGB-Diameter Curves')); +grid on; +box on; + +legend(traitp.tag,'Location','NorthWest'); \ No newline at end of file diff --git a/functions/plot_multicase_blmax.m b/functions/plot_multicase_blmax.m new file mode 100644 index 0000000..5e4549f --- /dev/null +++ b/functions/plot_multicase_blmax.m @@ -0,0 +1,35 @@ +function plot_multicase_blmax(dbh,blmax,traitp,fign) + +[ncases,~] = size(dbh); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + + +[nsets,ntimes] = size(blmax); + +usesets = 1:nsets; + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,9,4.5],'Color','w'); +subplot(1,2,1); +phan=plot(dbh(usesets,:)',blmax(usesets,:)','LineWidth',2); +xlim([0,20]); +xlabel('DBH [cm]'); +ylabel('Maximum Leaf Biomass [kg]'); +title(sprintf('Allometry Visual Tests\nIntercomparison of Leaf Biomas and Diameter Curves')); +grid on; +box on; +legend(traitp.tag,'Location','NorthWest'); + +subplot(1,2,2); +phan=loglog(dbh(usesets,:)',blmax(usesets,:)','LineWidth',2); +xlim([0,max(max(dbh(usesets,:)))]); +xlabel('DBH [cm]'); +ylabel('Maximum Leaf Biomass [kg]'); +title(sprintf('Allometry Visual Tests\nIntercomparison of Leaf Biomas and Diameter Curves')); +grid on; +box on; diff --git a/functions/plot_multicase_blmax_o_dbagdh.m b/functions/plot_multicase_blmax_o_dbagdh.m new file mode 100644 index 0000000..d04b2e9 --- /dev/null +++ b/functions/plot_multicase_blmax_o_dbagdh.m @@ -0,0 +1,27 @@ +function plot_multicase_blmax_o_dbagdh(dbh,blmax_o_dhdbag,traitp,fign) + +[ncases,~] = size(dbh); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + + +[nsets,ntimes] = size(dbh); + + +blmax_o_dhdbag(blmax_o_dhdbag>1e7 | blmax_o_dhdbag<1e-7)=NaN; + +usesets = 1:nsets; + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,7,6],'Color','w'); +phan=plot(dbh(usesets,:)',blmax_o_dhdbag(usesets,:)','LineWidth',2); +xlabel('DBH [cm]'); +ylabel('bl / dAGB/dH [m]'); +title(sprintf('Leaf Biomas per Relative Carbon Gain Rate per Height Rate ')); +grid on; +box on; +legend(traitp.tag,'Location','NorthWest'); diff --git a/functions/plot_multicase_blmaxdi.m b/functions/plot_multicase_blmaxdi.m new file mode 100644 index 0000000..29fcdf1 --- /dev/null +++ b/functions/plot_multicase_blmaxdi.m @@ -0,0 +1,20 @@ +function plot_multicase_blmaxdi(blmaxi,blmaxd,traitp,fign) + +[ncases,~] = size(blmaxi); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(blmaxi',blmaxd','LineWidth',2); +xlabel('blmax (direct) [kgC]'); +ylabel('blmax (integrated) [kgC]'); +title(sprintf('Allometry Visual Tests\nIntercomparison of Integrated and Direct Leaf Biomass')); +grid on; +box on; + +legend(traitp.tag,'Location','SouthEast'); \ No newline at end of file diff --git a/functions/plot_multicase_bsap.m b/functions/plot_multicase_bsap.m new file mode 100644 index 0000000..628e4e1 --- /dev/null +++ b/functions/plot_multicase_bsap.m @@ -0,0 +1,20 @@ +function plot_multicase_bsap(dbh,bsap,traitp,fign) + +[ncases,~] = size(bsap); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(dbh',bsap','LineWidth',2); +xlabel('dbh [cm]'); +ylabel('bsap [cm]'); +title(sprintf('Allometry Visual Tests\nSapwood Biomas [kgC] per diameter')); +grid on; +box on; + +legend(traitp.tag,'Location','SouthEast'); \ No newline at end of file diff --git a/functions/plot_multicase_bsapbag.m b/functions/plot_multicase_bsapbag.m new file mode 100644 index 0000000..faafc0e --- /dev/null +++ b/functions/plot_multicase_bsapbag.m @@ -0,0 +1,20 @@ +function plot_multicase_bsapbag(dbh,bsapbag,traitp,fign) + +[ncases,~] = size(dbh); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(dbh',bsapbag','LineWidth',2); +xlabel('DBH [cm]'); +ylabel('Fraction of AGB in Sapwood'); +title(sprintf('Allometry Visual Tests\nIntercomparison of Sapwood Fraction')); +grid on; +box on; + +legend(traitp.tag,'Location','SouthEast'); \ No newline at end of file diff --git a/functions/plot_multicase_bsapid.m b/functions/plot_multicase_bsapid.m new file mode 100644 index 0000000..7bf0ab0 --- /dev/null +++ b/functions/plot_multicase_bsapid.m @@ -0,0 +1,20 @@ +function plot_multicase_bsapid(bsapd,bsapi,traitp,fign) + +[ncases,~] = size(bsapd); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(bsapd,bsapi,'LineWidth',2); +xlabel('bsap (direct method) [kgC]'); +ylabel('bsap (integrated method) [kgC'); +title(sprintf('Sapwood Biomas [kgC] Integrator Check')); +grid on; +box on; + +legend(traitp.tag,'Location','SouthEast'); \ No newline at end of file diff --git a/functions/plot_multicase_dbhe.m b/functions/plot_multicase_dbhe.m new file mode 100644 index 0000000..6d6a362 --- /dev/null +++ b/functions/plot_multicase_dbhe.m @@ -0,0 +1,20 @@ +function plot_multicase_dbhe(dbh,dbhe,traitp,fign) + +[ncases,~] = size(dbh); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(dbh',dbhe','LineWidth',2); +xlabel('DBH [cm]'); +ylabel('Effective DBH [cm]'); +title(sprintf('Allometry Visual Tests\nIntercomparison of Effective Diameter and Diameter Curves')); +grid on; +box on; + +legend(traitp.tag,'Location','SouthEast'); \ No newline at end of file diff --git a/functions/plot_multicase_h.m b/functions/plot_multicase_h.m new file mode 100644 index 0000000..db0d89b --- /dev/null +++ b/functions/plot_multicase_h.m @@ -0,0 +1,34 @@ +function plot_multicase_h(dbh,h,traitp,fign) + +[ncases,~] = size(dbh); + +cmap = cbrewer('qual','Set2',ncases); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,6,4.5],'Color','w'); +phan=plot(dbh',h','LineWidth',2); +xlabel('DBH [cm]'); +ylabel('Height [m]'); +title(sprintf('Allometry Visual Tests\nIntercomparison of Diameter Height Curves')); +grid on; +box on; + +legend(traitp.tag,'Location','SouthEast'); + + + + + + + + + + + + + +end \ No newline at end of file diff --git a/functions/plot_singlecase_cfractions.m b/functions/plot_singlecase_cfractions.m new file mode 100644 index 0000000..a43bfb9 --- /dev/null +++ b/functions/plot_singlecase_cfractions.m @@ -0,0 +1,59 @@ +function plot_singlecase_cfractions(dbh,blmax,bfrmax,bcr,bsap,bdead,traitp,ic,fign) + + +[ndbh] = numel(blmax); + + +cmap = cbrewer('qual','Set2',4); +set(0,'defaultAxesColorOrder',cmap); +set(0,'defaultAxesFontSize',13); +set(0,'defaultAxesLineWidth',1); + + +btot = blmax+bfrmax+bdead+bsap; + +afrac = zeros(ndbh,4); +atot = zeros(ndbh,4); + +afrac(:,4) = bdead./btot; +afrac(:,3) = bsap./btot; +afrac(:,1) = blmax./btot; +afrac(:,2) = bfrmax./btot; +atot(:,4) = bdead; +atot(:,3) = bsap; +atot(:,1) = blmax; +atot(:,2) = bfrmax; +atags = {'blmax','bfrmax','bsap','bdead'}; + + +figure(fign); +clf; +set(fign,'Units','inches','Position',[fign,fign,7.5,4],'Color','w'); +subplot(1,2,1); +ahan=area(dbh,afrac); +for ihan=1:length(ahan) + set(ahan(ihan),'FaceColor',cmap(ihan,:)); +end +set(gca,'Xscale','log'); +xlabel('dbh [cm] (logscale)'); +ylabel('[-]'); +ylim([0,1]); +title(sprintf('Fraction of Total Carbon Allocated\nCase: %s',traitp.tag{ic})); +grid on; +box on; +legend(atags,'Location','SouthEast'); + +subplot(1,2,2); +syhan = semilogy(dbh,atot,'LineWidth',2); +for ihan=1:length(syhan) + set(syhan(ihan),'Color',cmap(ihan,:)); + if(strcmp(atags{ihan},'bfrmax')) + display('brow'); + set(syhan(ihan),'LineStyle','--'); + end +end +title(sprintf('Total Carbon Allocated\nCase: %s',traitp.tag{ic})); +grid on; +box on; +xlabel('dbh [cm]'); +ylabel('[kgC]'); diff --git a/gen_param_instance.m b/gen_param_instance.m new file mode 100644 index 0000000..a490657 --- /dev/null +++ b/gen_param_instance.m @@ -0,0 +1,49 @@ +function [pdat] = gen_param_instance + +% Manually Create a Group of Parameters + +pdat.tag = {'def-chave14-lecure', 'def-chave14-calvocarapa','def-poorter06-lescure','2par-pine', '2par-maple', 'def-SaldaObrien' , 'King-SaldaObrien'}; +pdat.h_mode = {'chave14', 'chave14', 'poorter06', '2par_pwr', '2par_pwr', 'obrien' , 'obrien'}; +pdat.blmax_mode = {'2par_pwr_asym', '2par_pwr_asym', '2par_pwr_asym', 'salda', 'salda', 'salda' , 'salda'}; +pdat.bfrmax_mode = {'constant', 'constant', 'constant', 'constant', 'constant', 'constant' , 'constant'}; +pdat.bag_mode = {'chave14', 'chave14', 'chave14', '2par_pwr', '2par_pwr', 'salda' , 'salda'}; +pdat.bcr_mode = {'constant', 'constant', 'constant', 'constant', 'constant', 'constant' , 'constant'}; +pdat.bsap_mode = {'constant', 'constant', 'constant', 'constant', 'constant', 'constant' , 'constant'}; +pdat.wood_density = [0.7, 0.7, 0.7, 0.39, 0.56, 0.7 , 0.7]; +pdat.latosa_int = [0.75, 0.75, 0.75, 0.75, 0.75, 0.75 , 0.75]; +pdat.latosa_slp = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0]; +pdat.c2b = [2.0, 2.0, 2.0, 2.0, 2.0, 2.0 , 2.0]; +pdat.eclim = [-0.15, 0.0, -999, -999, -999, -999 , -999]; +pdat.bl_min = [0.4, 0.4, 0.4, 0.4, 0.4, -999 , -999]; + +pdat.froot_leaf = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0 , 1.0]; +pdat.ag_biomass = [0.7, 0.7, 0.7, 0.7, 0.7, 0.7 , 0.7]; + +pdat.h_max = [35.0, 35.0, 61.0, 35.0, 35.0, 35.0 , 35.0]; +pdat.h_min = [0.5, 2.5, 0.5, 0.5, 0.5, 2.5 , 2.5]; + +pdat.slatop = [0.024, 0.024, 0.024, 0.008, 0.03, 0.024 , 0.024]; + +pdat.dbh_adult = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 10.0]; +pdat.dbh_sap = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0]; + +pdat.d2h1 = [0.893, 0.893, 61.7, exp(0.58119), exp(0.31414), 0.64, 0.64]; +pdat.d2h2 = [0.760, 0.760, -0.0352, 0.56063, 0.71813, 0.37, 0.37]; +pdat.d2h3 = [-0.0340, -0.0340, 0.694, -999, -999, -999, -999]; + +pdat.d2h1_sap = [-999, -999, -999, -999, -999, -999, 1.0729]; +pdat.d2h2_sap = [-999, -999, -999, -999, -999, -999, 1.4925]; + +pdat.d2bl1 = [0.00873, 0.012, 0.00873, 0.0419, 0.0419, 0.0419 , 0.0419]; +pdat.d2bl2 = [2.1360, 2.089, 2.1360, 1.56, 1.56, 1.56 , 1.56]; +pdat.d2bl3 = [-999, -999, -999, 0.55, 0.55, 0.55 , 0.55]; + +pdat.d2bl1_sap = [-999, -999, -999, -999, -999, -999, 0.0201]; +pdat.d2bl2_sap = [-999, -999, -999, -999, -999, -999, 3.1791]; + + +pdat.d2bag1 = [0.0673, 0.0673, 0.0673, exp(-2.6177), exp(-1.801), -999 , -999]; +pdat.d2bag2 = [0.976, 0.976, 0.976, 2.4638, 2.3852, -999 , -999]; + + +end \ No newline at end of file