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!