Track down error in R

[ng.hanh]

Hi everyone!
I am new to modelling in R and sometimes I have code for the model with too many variables and it is really killing me to come after the error to fix. I wondering if someone knows the best way to find bugs and fix. For example, I've been working on this codes for 2 weeks.

library(deSolve)
#rm(list=ls()) #clear the workspace

#function to create a list(vector) of parameters
get_parameters=function(){
  K=0.4 #part of C translocated to the rhizomes 
  epsilon2=1 #belowground N content limitation coeff. 
  RNquotmin=0.01 #minimum belowground N content
  RNquotmax=0.05 #maximum belowground N content
  delta1=0.6 #limitation coeff for leaf N uptake
  LVm=0.007 #max leaf N uptake rate (mol N(mol C)-1 per day)]
  KL=9.2 #half sat limit for leaf uptake (mmol N m-3)
  RECmax=0.039 #theoretical maximum recruitment rate (per day)
  Krec2=10 #half sat coeff for limitation by the belowground biomass (g DW m-2)
  SB0=0.05 #initial biomass of new shoot (mmol C)
  Orec=1.1 #recruitment increasing rate with T
  LMR20c=0.025 #max leaf mortality rate t 20c (per day)
  D=1.4 #water column depth (meters)
  min0=0.2 #N mineralization rate in water at 0 C (per day)
  Kmin=0.5 #half sat constant for mineralization (mg/m3)
  tt=27 #degrees celcius
  K1=0.4 #light extinction coeff due to water (m-1)
  O2water=5.35 #oxygen concentration in water column g/m-3
  O2sed=2.5 #oxygen concentration in sediment profile g/m-3
  Dsed=1 #sediment depth meters
  poro=0.8 #porosity of sediment 
  w=1.1 #photosynthetic quotient (molO2(molC)-1)
  Pmax0c=-0.0467 #theoretical max production at 0 C (gO2(mmolC degreesC)-1 per day)
  KLAI=0.000864 #LAI/LB ratio (m2(mmolC)-1)
  I=600 #measured value - PAR at sea surface (W m-2)
  K2=4.8 #light ext. coeff. due to epiphytes
  K3=0.6 #absorption coeff. linked to canopy optical and geometrical properties 
  Ikmin=35 #min saturation light intensity (W m-2)
  Ikmax=80 #max saturation light intensity (W m-2)
  d=1 #day in the year
  epsilon1=0.4 #leaf N content limitation coeff 
  Opmax=0.0265 #Production increasing rate with T (gO2(mmolC degreesC)-1 per day)
  OLM=1.1 #leaf mortality increasing rate with T
  LMRv=0.08 #leaf sloughing coeff due to wind
  Vvent= #wind speed (m s-1)
  K4=1.2 #wind effect attenuation with depth (m-1)
  EGRmax=0.34 #max epiphytes growth (per day)
  IK2=40 #sat light intensity (W m-2)
  Oe=1.1 #growth increasing rate with T
  Ke=2 #half saturation coeff for the N limitation (mmolN m-3)
  KLGR=0.1 #half sat coeff for the limitation by leaf growth rate (per day)
  Srec=0.9 #internal reclamation threshold
  Ecn=9 #C:N ratio for epiphytes (molC(molN)-1)
  OLR=0.000045 #leaf respiration increasing rate with T (gO2(mmolC degreeC)-1 per day)
  ORR=0.000014 #rhizomes/root respiration increasing rate with T (gO2(mmolC degreeC)-1 per day)
  LR0c=0.00059 #theoretical leaf respiration at 0c (gO2(mmolC)-1 per day)
  RR0c=0.000147 #theoretical rhizomes/root respiration at 0c (gO2(mmolC)-1 per day)
  tau=0.1 #N transfer speed btw leaves and rhizomes (molN(molC)-1 per day)
  RMR20c=0.025 #Max rhizomes/root mortality rate at 20c (per day)
  delta2=0.6 #limitation coeff for root N uptake
  Kbnit=2 #half sat constant for nit/denit (mg m-3)
  bnit0=0.1 #nit rate at 0c (d-1)
  denit0=0.3 #denit rate at 0c (d-1)
  KT=0.07 #T coeff (degreesC-1)
  Orec=1.1 #recruitment increasing rate with T
  LNquotmin=0.03
  LNquotmax=0.07
  RVm=0.0035 #max N root uptake rate (molN(molC)-1 per day)
  KR=104 #half sat coeff for root uptake (mmolN m-3)
  
  parnames=c(ls())
  npars=length(parnames)
  pars=numeric(npars)
  for (i in 1:npars){pars[i]=(get(parnames[i]))}
  names(pars)=parnames
  return(pars)
}

dydt = function(t,y,pars,EGR,EN,REC,LM,LGR, LABSNH4,LABSNO3,Rtrans,Ttrans,LNrec,min,EABSNH4,EABSNO3,bmin,RABSNH4,bnit,denit){
  with(as.list(pars),{
    #Aboveground pools
    EB=y[1] #epiphytes biomass (mmolC m-2)
    SD=y[2] #density of shoots (m-2) 
    LB=y[3] #aboveground biomas (mmol C m-2)
    LN=y[4] #aboveground N pool (mmol N m-2)
    Ndet=y[5] #organic particulate N pool - detritus
    NH41=y[6] #aboveground ammonium 
    NO31=y[7] #aboveground nitrate
    #Belowground pools
    RB=y[8] #belowground biomass (mmol C m-2)
    RN=y[9] #belowground N pool (mmol N m-2)
    NH42=y[10] #belowground N pool (mmol N m-2)
    NO32=y[11] #belowground nitrate
    
    #Aboveground pools
    #dEB functions
    f5tt=OLM^(tt-20) #Mortality limitation due to temp
    f6v=LMRv*(1/10*(Vvent))*exp(-K4*D) #leaf loss function due to wind henerated currents and waves (per day)
    LM=LMR20c*f5tt+f6v #leaf loss rate (per day)
    f1E=exp(-K2*(EB/LB))
    Qcan=I*f1E*exp(-K1*(D-0.3))
    f10l=tanh((Qcan/f1E)/IK2) #Light limitation function (dimensionless)
    f11tt=Oe^(tt-20) #T limitation function 
    f12N=NH41+NO31/NH41+NO31+Ke #DIN limitation function
    IK1=1/2*(Ikmax+Ikmin)+1/2*(Ikmax-Ikmin)*cos(2*pi*(1/365)*d-110)
    fz=function(z,I,f1E,K1,D,K3,IK1){
      Qcan<-I*f1E*exp(-K1*(D-0.3))
      f2l<-exp(-K3*z)
      fz<-tanh(Qcan*f2l/IK1)
    }
    LAI=KLAI*LB
    int1<-integrate(fz,lower=0,upper=LAI,I=I,f1E=f1E,K1=K1,D=D,K3=K3,IK1=IK1)
    LNquot=LN/LB #leaf N quota (mmol N(mmol C)-1)
    LNsat=(LNquot-LNquotmin)/(LNquotmax-LNquotmin) #sat level of leaf N quota
    f3LN=LNsat^epsilon1
    Pmax=Opmax*t-Pmax0c #max production rate (gO2(mmolC)-1 per day)
    Ptot=Pmax*int1[1]$value #total productivity 
    LR=OLR*tt+LR0c #leaf respiration (gO2(mmolC)-1 per day)
    RR=ORR*tt+RR0c #rhizomes/root respiration (gO2(mmolC)-1 per day)
    TG=(Ptot*f3LN-(LR+RR))*(LB*1000/32*w)
    if(TG<0){
      TG<-0
    } #total net growth rate(mmol C m-2/day)
    RNquot=RN/RB #belowground N quota (mmol N(mmol C)-1)
    RNsat=(RNquot-RNquotmin)/(RNquotmax-RNquotmin) #saturation level of N belowground quota
    f4RN=RNsat^epsilon2 #N limitation function for rhizomes/roots
    RGR=TG*f4RN*K #belowground growth rate 
    LGR=TG-RGR #leaf growth rate (mmol C m-2 per day)
    f13LG=1-((LGR/LB)/(LGR/LB)+KLGR) #leaf growth rate limitation function (dimensionless)
    EGR=EGRmax*f10l*f11tt*f12N*f13LG
    LM=LMR20c*f5tt*f6v #leaf loss rate per day 
    EM=LM
    dEB=(EGR-EM)*EB
    
    #dSD functions
    f7tt=Orec^(tt-12) #theoretical limitation function for recruitment (function for T between 5-22 C)
    #dRB functions
    f9RB=RB/(RB+Krec2) #belowground limitation for recruitment
    REC=RECmax*f7tt*f9RB #shoot recruitment rate (per day)
    RM=RMR20c*f5tt
    dRB=RGR-REC*SB0*SD-RM*RB
    dSD=(REC-LM)*SD
    
    #dLB functions
    dLB=LGR+REC*SB0*SD-LM*LB
    
    #dLN functions
    LABSNH4=LVm*(NH41/NH41+KL)*(1-LNsat)^delta1
    LABSNO3=LVm*(NO31/NO31+KL)*(1-LNsat)^delta1
    #Rtrans - C transfer rate from rhizomes to leaves (mmolN(mmolC)-1) per day
    #Ltrans - C transfer rate from leaves to rhizomes (mmolN(mmolC)-1) per day
    deltasat=RNsat-LNsat
    if(deltasat>0&LNsat<0.75){
      Ltrans<-0
      Rtrans<-tau*abs(deltasat)
    }else if(deltasat<0&RNsat<0.75){
      Ltrans<-tau*abs(deltasat)
      Rtrans<-0
    }else{
      Ltrans<-0
      Rtrans<-0
    }
    #LNrec - part of N reclaimed inside the roots/rhizomes
    LNrec=(1-(LNsat/Srec)^2)*RECmax
    if(LNsat>Srec){
      LNrec<-0
    }
    dLN=(LABSNH4+LABSNO3)*LB+REC*SB0*SD*(RN/RB)+(Rtrans-Ltrans)*LN-(1-LNrec)*LM*LN
    
    #dNdet functions
    ftt=exp(KT*tt)
    min=min0*ftt*(O2water/O2water+Kmin) #mineralization in water column
    dNdet=(1-LNrec)*LM*(LN/D)+LM*(EB/9*D)-min*Ndet
    
    #dNH41 functions
    EABSNH4=(EGR/Ecn)*(NH41/NH41+NO31) #ammonium uptake by epiphytes (mmolN(mmolC)-1 per day)
    dNH41=min*Ndet-LABSNH4*(LB/D)-EABSNH4*(EB/D)
    
    #dNO31 functions
    EABSNO3=(EGR/Ecn)*(NO31/NH41+NO31) #nitrate uptake by epiphytes (mmolN(mmolC)-1 per day)
    dNO31=-LABSNO3*(LB/D)-EABSNO3*(EB/D)
    
    #dNH42 functions
    bmin=min
    RABSNH4=RVm*(NH42/NH42+KR)*(1-RNsat)^delta2 #ammonium uptake by the roots (mmolN(mmolC)-1 per day)
    dNH42=bmin*(1-poro/poro)*Ndet-(RABSNH4/Dsed*poro)*RB
    
    #dRN functions
    RNrec=(1-(RNsat/Srec)^2)*RECmax
    if(RNsat>Srec){
      RNrec<-0
    }
    dRN=RABSNH4*RB-REC*SB0*SD*(RN/RB)+(Ltrans-Rtrans)*RN-(1-RNrec)*RM*RN
    
    #dNO32 functions
    bnit=bnit0*ftt*(O2sed/O2sed+Kbnit)
    denit=denit0*ftt*(1-(O2sed/O2sed+Kbnit))
    dNO32=bnit*NH42-denit*NO32
    
    return(list(c(dEB,dSD,dLB,dLN,dNdet,dNH41,dNO31,dRB,dNH42,dRN,dNO32)))
  })
}

pars = get_parameters()

#get initial values of state variables - from equilibrium after running model to equilibrium

  EB <- 1
  SD <- 1
  LB <- 1
  LN <- 1
  Ndet <- 1
  NH41 <- 1
  NO31 <- 1
  RB <- 1
  NH42 <- 1
  RN <- 1
  NO32 <-1

yini = c(EB,SD,LB,LN,Ndet,NH41,NO31,RB,NH42,RN,NO32)
t=seq(0,100,0.01)#time (days); time step 0.01 days

#modify rates of input
out=lsoda(yini,times=t,func=dydt,parms=pars)

Thank you so much!

[mara]

Hi,

It's quite hard to read your code right now, as it's not formatted.

Ideally, you'd please turn this into a self-contained reprex (short for minimal reproducible example). However, if, for some reason, you cannot, check out the code-formatting FAQ:

FAQ: How to format your code

Formatting Code: the Basics All code or console output you include in your posts should be formatted properly. Luckily, this is very easy to do! Just use the code formatting button at the top of the post editing box: [image] Select some code Click the </> button! The code formatting button automatically adds special Markdown formatting symbols to whatever you selected. If you check the preview pane to the right of your post, you'll see that the forum software now outputs nicely formatted code . [04] The code formatting button is pretty smart! It can tell the difference between snippets of code that are part of a sentence and blocks of code that stand alone. …

Right now the best way to install reprex is:

# install.packages("devtools")
devtools::install_github("tidyverse/reprex")

If you've never heard of a reprex before, you might want to start by reading the tidyverse.org help page. The reprex dos and don'ts are also useful.

For pointers specific to the community site, check out the reprex FAQ, linked to below.

FAQ: What's a reproducible example (`reprex`) and how do I create one?

Why reprex? Getting unstuck is hard. Your first step here is usually to create a reprex, or reproducible example. The goal of a reprex is to package your code, and information about your problem so that others can run it and feel your pain. Then, hopefully, folks can more easily provide a solution. What's in a Reproducible Example? Parts of a reproducible example: background information - Describe what you are trying to do. What have you already done? complete set up - include any library() calls and data to reproduce your issue. data for a reprex: Here's a discussion on setting up data for a reprex make it run - include the minimal code required to reproduce your error on the data…