daily pastebin goal
83%
SHARE
TWEET

Untitled

a guest Oct 12th, 2018 52 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. model {
  2.   for(i in 2:nday){
  3.     dayl = 0
  4.     pi = 2*arcsin(1)
  5.     Tr = .5*exp(p10*.5*(mint[i] + maxt[i]))
  6.     lai = max(.1, Cf/LMA)
  7.     trange = (maxt[i] - mint[i])/2
  8.     gs = pow(2, 0.789798) / (0.37836)
  9.     pp = lai*Nit / gs *p11*exp(0.011136*maxt[i])
  10.     qq = -204.6453
  11.     Caqqpp.sq = pow((Ca[i] + qq - pp), 2)
  12.     ci = .5*(Ca[i] + qq - pp + pow( (Caqqpp.sq) - 4*(Ca[i]*qq - pp*4.22273), 1/2 ))
  13.    
  14.     e0 = 7.19298*pow(lai, 2) / (pow(lai, 2) + 2.1001)
  15.    
  16.     dec = -23.4*cos(360*(yearday+10)/365 *pi/180)*pi / 180
  17.     mult = tan(dec)*tan(lat)
  18.    
  19.     dayl = ifelse(abs(mult) < 1, 24*arccos(-mult)/pi, ifelse(mult >= 1, 24, 0) )
  20.    
  21.     cps = e0*rad[i]*gs*(Ca[i] - ci) / (e0*rad[i] + gs*(Ca[i] - ci))
  22.     G[i] = cps*(0.0156935*dayl + 0.0453194)
  23.    
  24.     Cf.pred[i] = (1-p5)*Cf[i-1] + G[i]*(1 - p2)*p3
  25.     Cw.pred[i] = (1 - p6)*Cw[i-1] + G[i]*(1-p2)*(1-p3)*(1-p4)
  26.     Cr.pred[i] = (1 - p7)*Cr[i-1] + G*(1-p2)*(1-p3)*p4
  27.     Clit.pred[i] = Cf[i-1]*p5 + Cr[i-1]*p7 + (1 - Tr*(p8+p1))*Clit[i-1]
  28.     Csom.pred[i] = p6*Cw[i-1] + Tr*p1*Clit[i-1] + (1 - Tr*p9)*Csom[i-1]
  29.    
  30.     Cf[i] ~ dnorm(Cf.pred[i], sd_cf)
  31.     Cw[i] ~ dnorm(Cw.pred[i], sd_cw)
  32.     Cr[i] ~ dnorm(Cr.pred[i], sd_cr)
  33.     Clit[i] ~ dnorm(Clit.pred[i], sd_clit)
  34.     Csom[i] ~ dnorm(Csom.pred[i], sd_csom)
  35.   }
  36.  
  37.   for(i in 1:nday){
  38.     Cf.obs[i] ~ dnorm(Cf[i], sd_cf_obs)
  39.     Cw.obs[i] ~ dnorm(Cw[i], sd_cw_obs)
  40.     Cr.obs[i] ~ dnorm(Cr[i], sd_cr_obs)
  41.     Clit.obs[i] ~ dnorm(Clit[i], sd_clit_obs)
  42.     Csom.obs[i] ~ dnorm(Csom[i], sd_csom_obs)
  43.   }
  44.  
  45.   ## Initial conditions
  46.   Cf[1] ~ dnorm(100, 2)
  47.   Cw[1] ~ dnorm(9200, 25)
  48.   Cr[1] ~ dnorm(100, 2)
  49.   Clit[1] ~ dnorm(20, 2)
  50.   Csom[1] ~ dnorm(11000, 25)
  51.  
  52.   ## Priors on parameters
  53.   p1 ~ dunif((10^(-6)), .01)
  54.   p2 ~ dunif(.2, .7)
  55.   p3 ~ dunif(.01, .5)
  56.   p4 ~ dunif(.01, .5)
  57.   p5 ~ dunif((10^(-4)), .1)
  58.   p6 ~ dunif((10^(-6)), .01)
  59.   p7 ~ dunif((10^(-4)), .01)
  60.   p8 ~ dunif((10^(-5)), .1)
  61.   p9 ~ dunif((10^(-6)), .01)
  62.   p10 ~ dunif(.05, .2)
  63.   p11 ~ dunif(5, 20)
  64.  
  65.   ## Priors on process errors
  66.   sd_cf ~ dnorm(0, 100)T(0,)
  67.   sd_cw ~ dnorm(0, 100)T(0,)
  68.   sd_cr ~ dnorm(0, 100)T(0,)
  69.   sd_clit ~ dnorm(0, 100)T(0,)
  70.   sd_csom ~ dnorm(0, 100)T(0,)
  71.  
  72.   sd_cf_obs ~ dnorm(0, 100)T(0,)
  73.   sd_cw_obs ~ dnorm(0, 100)T(0,)
  74.   sd_cr_obs ~ dnorm(0, 100)T(0,)
  75.   sd_clit_obs ~ dnorm(0, 100)T(0,)
  76.   sd_csom_obs ~ dnorm(0, 100)T(0,)
  77. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top