hi sir, actually I start receiving a new error (Error in eval(e, x, parent.frame()) : object 'DOYS' not found)
in the following code.
Lintul_1 <-function(Time, State, Pars, WDATA){
with(as.list(c(State, Pars)), {
#Get the weather data for the day
WDATA <- subset(WDATA, DOYS == floor(Time))
if(nrow(WDATA) != 1){
print(paste0("ERROR in weather data file: no or more than one record for day nr: ",Time))
print(WDATA)
stop()
}
DTEFF <- max(0, WDATA[,"DAVTMP"] - TBASE) # effective daily temperature (for crop development a threshold temperature (TBASE) needs to be exceeded)
RPAR <- FPAR * WDATA[,"DTR"] # PAR MJ m-2 d-1
#determine rates when crop is still growing
if(TSUM < FINTSUM){
# Once the emergence date is reached and enough water is available the crop emerges (1),
# Once the crop is established is does not disappear again (2)
emerg1 <- ifelse((Time - DOYEM) >= 0, 1, 0)
emerg2 <- ifelse(LAI > 0, 1, 0)
# Emergence of the crop is used to calculate the accumulated temperature sum.
EMERG <- max(emerg1,emerg2)
RTSUM <- DTEFF * EMERG
# Light interception (MJ m-2 d-1) and total crop growth rate (g m-2 d-1)
PARINT <- RPAR * (1 - exp(-K * LAI))
GTOTAL <- LUE * PARINT
# Relative death rate (d-1) due to aging
RDRDV <- ifelse((TSUM-TSUMAN) < 0, 0, approx(RDRT.TEMP,RDRT.RDR,WDATA[,"DAVTMP"])$y)
# Relative death rate (d-1) due to self shading
RDRSH <- pmax(0, pmin( RDRSHM * (LAI-LAICR) / LAICR, RDRSHM))
# Effective relative death rate (1; d-1) and the resulting decrease in LAI (2; m2 m-2 d-1) and leaf weight (3; g m-2 d-1)
RDR <- max(RDRDV, RDRSH) # (1)
DLAI <- LAI * RDR # (2)
RWLVD <- WLVG * RDR # (3)
# Allocation to roots (2), leaves (4), stems (5) and storage organs (6)
FRT <- approx(x = FRTTB.TSUM,y = FRTTB.FRT,xout = TSUM)$y # (2)
FLV <- approx(x = FLVTB.TSUM,y = FLVTB.FLV,xout = TSUM)$y # (4)
FST <- approx(x = FSTTB.TSUM,y = FSTTB.FST,xout = TSUM)$y # (5)
FSO <- approx(x = FSOTB.TSUM,y = FSOTB.FSO,xout = TSUM)$y # (6)
# Change in biomass (g m-2 d-1) for leaves (1), green leaves (2), stems (3), storage organs (4) and roots (5)
RWLV <- GTOTAL * FLV # (1)
RWLVG <- GTOTAL * FLV - RWLVD # (2)
RWST <- GTOTAL * FST # (3)
RWSO <- GTOTAL * FSO # (4)
RWRT <- GTOTAL * FRT # (5)
GLAI <- gla(TIME = Time, DTEFF = DTEFF, TSUM = TSUM, LAI = LAI, RWLV = RWLV, PARS = Pars)
# Change in LAI (m2 m-2 d-1) due to new growth of leaves
RLAI <- GLAI - DLAI
RATES <- c(RTSUM = RTSUM, RPAR = RPAR, RLAI = RLAI,
RWLV = RWLV, RWLVG = RWLVG, RWLVD = RWLVD,
RWST = RWST, RWSO = RWSO, RWRT = RWRT)
AUX <- c(WBIOMTHA = (WSO + WLV + WST + WRT) * 0.01,
WSOTHA = WSO * 0.01,
HI = WSO /(WSO + WLV + WST + WRT))
}
else{
#all plant related rates are set to 0
RATES <- c(RTSUM = 0, RPAR = 0, RLAI = 0,
RWLV = 0, RWLVG = 0, RWLVD = 0,
RWST = 0, RWSO = 0, RWRT = 0)
AUX <- c(WBIOMTHA = (WSO + WLV + WST + WRT) * 0.01,
WSOTHA = WSO * 0.01,
HI = ifelse((WSO + WLV + WST + WRT)==0, 0, WSO /(WSO + WLV + WST + WRT)))
}
return(list(RATES,c(AUX, RATES) ))
})
}
This function reads in CABO format weather files. If end time is >365 or >366 days, two years are combined.