;recreate farquar solver for faster solution? let's see


FUNCTION ed_calc_mech,ppm,vm,t,rh,par,c4=c4

  COMMON edmech, mech

;setup common struct
  mech = {lambda: 4.5e-2, theta: 38.4, D0: 0.01, kappa: -0.5, mol2gC: 12.01, mol2gh2o: 18.0153, M: 8.0, b: 1.0e4, alpha: 0.08, lgamma: 0.02, ca: 360.*1e-6, vm: 25., t: 12., ea: 0.5, L: 2000., Rn: 1000., c4: keyword_set(c4), ec:0., eo:0., ac:0., ao:0., tl:0., k1:0., k2:0., gamma:0., vmt:0., voff:35.}
  IF keyword_set(c4) THEN BEGIN
    mech.M = 4.
    mech.b = 0.01 * 1e6
    mech.alpha = 0.06
    mech.lgamma = 0.04
    mech.voff = 100.
  ENDIF 

;output is:
;0 1 2 3 4 5 = vm, ao, eo, tl, c1, cs

;input vars
  l = ((par>0.<3000.)/4.6)*4.72*(1.-exp(mech.kappa))
  rn = ((par>0.<3000.)/2.07)*(1.-exp(mech.kappa))
  ea =  exp(-5415.0/((t>(-12.)<60.)+273.15))*2.5414e6*((rh>10.<95.)/100.)
  ca = (ppm>100.<1000.) * 1.e-6
  
;initial values
  el = 2.5414e6 * exp(-5415.0/((t>(-12.)<60.)+273.15))
  ds = el-ea
  tl = ((rn - (mech.lambda*mech.b*ds))/mech.theta)+(t>(-12.)<60.)
  ci = ca * 0.7
  mech.vm = vm

  output = make_array(4,n_elements(ppm),/float)
  itmax=250000.
  FOR i = 0l,n_elements(ppm)-1l DO BEGIN 
    IF i MOD 10000l EQ 0l THEN print,'on elem ',i
    mech.T = t[i]
    mech.L = l[i]
    mech.Rn = rn[i]
    mech.ea = ea[i]
    mech.ca = ca[i]
    tlcase = tl[i]
    tlres = newton(tlcase,'ed_farqfunc_tl',/double,itmax=2500.)
    mech.tl = tlres
;    dtl = (1./288.2) - (1./(mech.tl+273.15))
    dtl = (1./298.2) - (1./(mech.tl+273.15))
;    mech.vmt = mech.vm * exp(3000.*dtl)/(1.+exp(0.4*(5.-mech.tl)))/(1.+exp(0.4*(mech.tl-mech.voff)))
 mech.vmt = mech.vm * exp(3000.*dtl)/(1.+exp(0.4*(10.-mech.tl)))/(1.+exp(0.4*(mech.tl-mech.voff)))
    mech.ec = mech.b * ((2.5414e6 * exp(-5415.0/(mech.tl+273.15)))-ea[i])
    mech.ac = (-1.0) * mech.vmt * mech.lgamma
    IF mech.L LT 10 THEN BEGIN
      mech.eo = mech.ec
      mech.ao = mech.ac
    ENDIF ELSE BEGIN 
      mech.gamma = 2.12e-5 * exp(5000.*dtl)
      mech.k1 = (1.5e-4) * exp(6000.*dtl)
      mech.k2 = 0.836 * exp(-1400.*dtl)
      IF keyword_set(c4) THEN BEGIN
        ao_init = ((min([mech.alpha*mech.l,mech.vmt,18000.*mech.vmt*ci[i]]))>0.) - (mech.vmt*mech.lgamma)
      ENDIF ELSE BEGIN 
        ao_init = ((min([mech.alpha * mech.l * ((ci[i]-mech.gamma)/(ci[i]+mech.gamma)),mech.vmt * (ci[i]-mech.gamma) / (ci[i] + (mech.k1*(1.+mech.k2)))]))>0.) - (mech.vmt * mech.lgamma)
      ENDELSE 
      cs = (1.6*ao_init)/(mech.ca-ci[i])
      eo_init = cs * ((2.5414e6 * exp(-5415.0/(mech.tl+273.15)))-ea[i])
      ac_init = (mech.b/1.6)*(mech.ca-ci[i])
 ;   clcase = [ac_init>(-5.)<(50.)]
 ;   clres = newton(clcase,'ed_farqfunc_cl',/double)
;    clres = broyden(clcase,'ed_farqfunc_cl')
 ;   mech.ac = clres
      opcase = [ao_init>(-5.)<(50.),eo_init>(1.)<(10000.)] ;[ao_init,ec_init]
;      opcase = newton(opcase,'ed_farqfunc',itmax=itmax)
      opres = newton(opcase,'ed_farqfunc',/double,itmax=itmax,tolf=1.e-3,tolx=1.e-6,tolmin=1.e-5)
;    opres = broyden(opcase,'ed_farqfunc')
      mech.ao = opres[0]
      mech.eo = opres[1]
;    print,mech.ao,mech.ac,mech.eo,mech.ec
    ENDELSE 
    output[*,i] = [mech.ao,mech.ac,mech.eo,mech.ec]
  ENDFOR 
  return,output
END


FUNCTION ed_farqfunc_tl,x
  COMMON edmech
  tl = x[0]
  el = 2.5414e6 * exp(-5415.0/(tl+273.15))
  ds = el - mech.ea
  eo = mech.b * ds
  return,[(mech.lambda*eo)+(mech.theta*(tl-mech.t))-mech.rn]
END

FUNCTION ed_farqfunc_cl,x
  COMMON edmech
  ac = x[0]
  ci = mech.ca-((1.6*ac)/mech.b)
  IF mech.c4 NE 0 THEN BEGIN 
    test = (((min([mech.alpha*mech.l,mech.vmt,18000.*mech.vmt*ci]))>0.)-(mech.vmt*mech.lgamma))-ac
  ENDIF ELSE BEGIN 
    test = (((min([mech.alpha * mech.l * ((ci-mech.gamma)/(ci+mech.gamma)),mech.vmt * (ci-mech.gamma) / (ci + (mech.k1*(1.+mech.k2)))]))>0.) - (mech.vmt * mech.lgamma))-ac
  ENDELSE 
  return,[test]
END 

FUNCTION ed_farqfunc,x
  COMMON edmech

;x = ao, eo
  ao = x[0] ;< (-10.) < 100.
  eo = x[1] ;> 10. < 20000.
 
  tl = ((mech.Rn - (mech.lambda * eo))/mech.theta)+mech.t
  tl>=-40.
  tl<=200.
  el = 2.5414e6 * exp(-5415.0/(tl+273.15))
  ds = el - mech.ea
;  dtl = (1./288.2) - (1./(Tl+273.15))
  dtl = (1./298.2) - (1./(Tl+273.15))
  gamma = 2.12e-5 * exp(5000.*dtl)
  k1 = (1.5e-4) * exp(6000.*dtl)
  k2 = 0.836 * exp(-1400.*dtl)
  vm = mech.vm * exp(3000.*dtl)/(1.+exp(0.4*(10.-tl)))/(1.+exp(0.4*(tl-mech.voff)))
;  vm = mech.vm * exp(3000.*dtl)/(1.+exp(0.4*(5.-tl)))/(1.+exp(0.4*(tl-mech.voff)))
  IF ~finite(vm) THEN stop
  IF ds EQ 0 THEN BEGIN
    cs = 0.
    ci = mech.ca
  ENDIF ELSE BEGIN 
    cs = Eo/ds
    ci = mech.ca - (1.6*ao/cs)
  ENDELSE 
  
  test1 =((mech.M * ao)/((ci-gamma)*(1.+(ds/mech.d0)))) + mech.b - cs
  IF mech.c4 NE 0 THEN BEGIN 
    test2 = ((min([mech.alpha * mech.l,vm,18000.*vm*ci])>0.)-(vm*mech.lgamma))-ao
  ENDIF ELSE BEGIN 
    test2 = (((min([mech.alpha * mech.l * ((ci-gamma)/(ci+gamma)),vm * (ci-gamma) / (ci + (k1*(1.+k2)))]))>0.) - (vm * mech.lgamma))-ao
  ENDELSE 

  return,[test1,test2]
END