function [o18e,o18l,o18s,o18c,o18cc,mcerr,o18ccp]=imc(pr,ta,rh,dfs,dfb,efc,blc,c10,c11,pl,feo,fwxm,c1,c2,c3,c4,c5,c6,c7,c8,c9)
% function [o18e,o18l,o18s,o18c,o18cc,mcerr]=imc(pr,ta,rh,dfs,dfb,efc,blc,c10,c11,pl,feo,fwxm,c1,c2,c3,c4,c5,c6,c7,c8,c9)
% Barbour-modified Roden model
% for use in constructing synthetic 
% oxygen isotope series 
% for comparison with collected data
% after Barbour et al. (2004) 
% (hereafter referred to as BRFE2004)
% See end of script for full references cited
%
% 
%
% default values for BRFE2004 model parameters are:
%
% dfs=32; % diffusive fractionation  through stomata (per mil)
% dfb=21; % diffusive fractionation  through boundary layer (per mil)
% efc=27; % equilibrium fractionation between C=O and water (per mil)
% blc=2100; % boundary layer conductance (mmol m-2 s-1)
% pl=0.018; % effective length for the Peclet effect (m)
% feo=0.42; % proportion of exchangeable oxygen in cellulose
% fwxm= 1.0; % proportion of xylem water in meristem
%
% parameters for timeseries model extension are:
%
% ta = air temperature(degC)
% pr = precipitation(mm/mo)
% rh = relative humidity(%)
% c1 = y-intercept of o18sw vs. precipitation amount relationship
% c2 = slope of o18sw vs. precipitation amount relationship
% c3 = y-intercept of leaf T vs. air T relationship
% c4 = slope of leaf T vs. air T relationship
% c5 = height(Km) at which precip reaching ground is condensing
% c6 = moist adiabatic lapse rate (K/Km)
% c7 = mean observed o18 of cellulose (for baseline correction)
% c8 = fraction RMS variance due to environmental forcing
% c9 = percent uncertainty in model parameters for Monte Carlo
%            error estimation
% c10 = first constant for stomatal conductance function of air vapor pressure deficit
%            (mb mmol m-2 s-1 )
% c11 = second constant for stomatal conductance function of air vapor pressure deficit
%            (mmol m-2 s-1)
%
% default values of timeseries model extension parameters are:
% 
% c1 = 0.74 permil; % [Lachniet and Patterson, 2004 (submitted)]
% c2 = -0.0285 permil/(mm/mo); % [Lachniet and Patterson, 2004]
% c3 = 10 degC; % [Linacre, 1964]
% c4 = 23/33 degC/degC; % [Linacre, 1964]
% c5 = 3 Km; % [after Lawrence and Gedzelman (2003); Lawrence et al. (2004)]
% c6 = 6 K/Km; % [e.g. Wallace and Hobbs, 1977]
% c7 = 27 permil; % [reasonable values are 24-35 permil]
% c8 = 0.45; % free parameter chosen to match variance in simulated
%              vs. observed data
% c9 = 20; % percent uncertainty in model parameters for Monte Carlo
%            error estimation

%error(nargchk(1,21,nargin));
if nargin <= 21,
  dfs=32; 
  dfb=21; 
  efc=27;
  blc=2100;
  pl=0.018;
  feo=0.42;
  fwxm= 1.0;
  c1 = 0.74;
  c2 = -0.0285;
  c3 = 10;
  c4 = 23/33;
  c5 = 3;
  c6 = 6;
  c7 = 27.17;
  c8 = 0.45;
  c9 = 20;
  c10= 4290;
  c11=-43;
end                                                          

  c=[dfs dfb efc blc c10 c11 pl feo fwxm c1 c2 c3 c4 c5 c6 c7 c8 c9];

% simulation with selected parameter values:
  [o18e,o18l,o18s,o18c,o18cc]=ims(pr,ta,rh,dfs,dfb,efc,blc,c10,c11,pl,feo,fwxm,c1,c2,c3,c4,c5,c6,c7,c8);

% apply Monte Carlo parameter error estimation to get error bars:
  o18ep=zeros(1,length(pr));
  o18lp=zeros(1,length(pr));
  o18sp=zeros(1,length(pr));
  o18cp=zeros(1,length(pr));
  o18ccp=zeros(1,length(pr));
  a=0;
  v=[(100-c(18))/100:2*c(18)/(100*(length(c)-1)):(100+c(18))/100];
  cm=zeros(1,length(v));
  mcerr=zeros(length(pr),2*5);

  if c(18)==0,
       cm=c;
       [o18emc,o18lmc,o18smc,o18cmc,o18ccmc]=ims(pr,ta,rh,cm(1),cm(2),cm(3),cm(4),cm(5),cm(6),cm(7),cm(8),cm(9),cm(10),cm(11),cm(12),cm(13),cm(14),cm(15),c(16),cm(17));
       o18ep=o18emc';
       o18lp=o18lmc';
       o18sp=o18smc';
       o18cp=o18cmc';
       o18ccp=o18ccmc';
       cmp=cm;
  elseif c(18) ~=0,
    for i=1:1000
      m=randperm(length(v));
      cm=c.*v(m);
      [o18emc,o18lmc,o18smc,o18cmc,o18ccmc]=ims(pr,ta,rh,cm(1),cm(2),cm(3),cm(4),cm(5),cm(6),cm(7),cm(8),cm(9),cm(10),cm(11),cm(12),cm(13),cm(14),cm(15),c(16),cm(17));
      o18ep(i,:)=o18emc';
      o18lp(i,:)=o18lmc';
      o18sp(i,:)=o18smc';
      o18cp(i,:)=o18cmc';
      o18ccp(i,:)=o18ccmc';
      cmp(i,:)=cm;
    end
  end
    if c(18)~=0,                                                                                                              
  o18eps=sort(o18ep);
  o18lps=sort(o18lp);
  o18sps=sort(o18sp);
  o18cps=sort(o18cp);
  o18ccps=sort(o18ccp);

  mcerr=[o18eps(50,:)' o18eps(950,:)' o18lps(50,:)' o18lps(950,:)' o18sps(50,:)' o18sps(950,:)' o18cps(50,:)' o18cps(950,:)' o18ccps(50,:)' o18ccps(950,:)'];
end

  return