#' Finds solar radiation reaching the ground surface
#'
#' @param lat Latitude in degrees
#' @param hr Hour of the day, 0-24, decimal format
#' @param month Name of the month, capitalised (e.g. "January")
#' @param aspect Aspect in degrees
#' @param slope Slope in degrees
#' @param LAI Leaf Area Index
#' @param cloud Proportion cloud cover (0-1)
#'
#' @export
insol <- function(lat = -32, hr = 12, month = "December", aspect = 0, slope = 0, LAI = 3, cloud = 0){
radians <- function(deg) {(deg * pi) / (180)}
hAngle <- ifelse(((15*(12-hr))<0), (360+(15*(12-hr))),(15*(12-hr)))
decl <- dplyr::case_when(month == "March" ~ 0,
month == "April" ~ -11.75,
month == "May" ~ -20.35,
month == "June" ~ -23.50,
month == "July" ~ -20.35,
month == "August" ~ -11.75,
month == "September" ~ 0,
month == "October" ~ 11.75,
month == "November" ~ 20.35,
month == "December" ~ 23.50,
month == "January" ~ 20.35,
TRUE ~ 11.75)
terrainE <- (radians(cos(radians(hAngle-aspect))*slope))
zenith <- acos(sin(radians(abs(lat)))*sin(radians(decl))+cos(radians(abs(lat)))*cos(radians(decl))*cos(radians(hAngle)))-terrainE
shade <- 1-exp(-(2/(pi*tan(max(0.01,1.570796-zenith))))*LAI)
out <- max(0,((1-shade)*(1-cloud)*1000*cos(zenith)))
return(out)
}
# Finds saturation specific humidity from Stull (1988) Eqn 7.5.2c and d
QSat <- function(tAir = 293.46, pAir = 101180){
Q <- 0.622 * (611.2 * exp(17.67 * (tAir - 273.15) / (tAir - 29.66))) / pAir
return(Q)
}
# Finds water vapour flux
# Calculates RH at the litter surface from inversion of Nelson (1984) EMC equation
Ema <- function(m, nelsonA = 5.2, nelsonB = -19, conLitter = 0.0006, tAir = 293.46, vAir = 3.52, pAir = 101180, insolation = 0, rhAir = 58){
SpecHum <- rhAir/100 *QSat(tAir, pAir) # Specific humidity
tLeaf <- tAir + (8.15 - 2.25 * exp(-0.6 * vAir) - 0.0312 * tAir + (0.021 + (-0.04 + 0.0006 * tAir - 0.00000125 * tAir ^ 2) * exp(-0.6 * vAir)) * insolation) # temperature of the leaf
fA <- -((-nelsonA / nelsonB) - 0.001 * (tLeaf - 293.15)) * nelsonB
rhLeaf <- exp((-18.0153 / (1.9872 * tLeaf)) * exp(m * nelsonB + fA))
out <- 1.14 * conLitter * (rhLeaf * QSat(tLeaf, pAir) - SpecHum) # Water vapour flux from litter
return(out)}
#' Finds the moisture content of dead leaf litter at a given time-step
#'
#' Function uses the 'Single differential equation model' of Matthews et al (2010)
#'
#' @param m The starting moisture (proportion ODW)
#' @param nelsonA Constant from Nelson FMC model
#' @param nelsonB Constant from Nelson FMC model
#' @param tAir Atmospheric temperature (K)
#' @param vAir Wind speed (m/s)
#' @param pAir Atmospheric pressure (Pa)
#' @param rhAir Relative humidity (Percent)
#' @param sigma Surface area: volume ratio of litter particles (m2/m3)
#' @param rhoLitter Density of leaves (kg/m3)
#' @param conLitter Litter conductance (m/s). Using value from Matthews spreadsheet rather than paper
#' @param dt Seconds per time step
#' @param EPS Sensitivity of iteration
#' @param insolation Solar energy at the soil surface (kW/m2)
#'
#' @export
simplefmc <- function(m, tAir = 293.46, vAir = 3.52, pAir = 101180, rhAir = 58, dt = 3600, insolation = 150, rhoLitter = 550,
nelsonA = 5.2, nelsonB = -19, conLitter = 0.0006, sigma = 3000, EPS = 0.01)
{
SpecHum <- rhAir/100 *QSat(tAir, pAir) # Specific humidity
xStart <- m
xEnd <- m
Ftest <- 0 - sigma * Ema(m, nelsonA, nelsonB, conLitter, tAir, vAir, pAir, insolation, rhAir)
while(abs(Ftest) > 0.0001){
Ftest <- rhoLitter * (xStart - xEnd) / dt
Ftest <- Ftest - sigma * Ema(m, nelsonA, nelsonB, conLitter, tAir, vAir, pAir, insolation, rhAir)
if(abs(Ftest) > 0.0001){
temp <- xEnd
h <- EPS * abs(temp)
if(h == 0) {h <- EPS}
xEnd <- temp + h
m <- xEnd
h <- xEnd - temp
F2 <- rhoLitter * (xStart - xEnd) / dt - sigma * Ema(m, nelsonA, nelsonB, conLitter, tAir, vAir, pAir, insolation, SpecHum)
xEnd <- temp
m <- xEnd
Ftest <- -Ftest * h / (F2 - Ftest)
xEnd <- xEnd + Ftest
m <- xEnd
}
}
fmc <- m
}
#' Helper function for frameWeather
#'
#' @param clim
#'
tidyWeather <- function(clim) {
# Fill empty directions
if (clim$dPM[1] == 0 || is.null(clim$dPM[1])) {
fill <- "N"
} else {
fill <- clim$dPM[1]}
if (clim$dAM[1] == 0 || is.null(clim$dAM[1])) {
clim$dAM[1] <- fill
}
if (clim$dPM[1] == 0 || is.null(clim$dPM[1])) {
clim$dPM[1] <- clim$dAM[1]
}
for (r in 2:nrow(clim)) {
if (clim$dPM[r] == 0 || is.null(clim$dPM[r])) {
fill <- "N"
} else {
fill <- clim$dPM[r]}
if (clim$dAM[r] == 0 || is.null(clim$dAM[r])) {
clim$dAM[r] <- fill
}
if (clim$dPM[r] == 0 || is.null(clim$dPM[r])) {
clim$dPM[r] <- clim$dAM[r]
}
}
# Set initial values
clim <- clim %>%
mutate(tAM = tAM + 273.15,
tPM = tPM + 273.15,
tMin = tMin + 273.15,
tMax = tMax + 273.15,
wAM = wAM/3.6,
wPM = wPM/3.6,
cAM = cAM/8,
cPM = cPM/8,
MSLP = MSLP*100,
# Set directions to degrees
dAM = case_when(dAM == "N" ~ 0,
dAM == "NNE" ~ 22,
dAM == "NE" ~ 45,
dAM == "ENE" ~ 68,
dAM == "E" ~ 90,
dAM == "ESE" ~ 112,
dAM == "SE" ~ 135,
dAM == "SSE" ~ 158,
dAM == "S" ~ 180,
dAM == "SSW" ~ 202,
dAM == "SW" ~ 225,
dAM == "WSW" ~ 248,
dAM == "W" ~ 270,
dAM == "WNW" ~ 292,
dAM == "NW" ~ 315,
dAM == "NNW" ~ 338),
dPM = case_when(dPM == "N" ~ 0,
dPM == "NNE" ~ 22,
dPM == "NE" ~ 45,
dPM == "ENE" ~ 68,
dPM == "E" ~ 90,
dPM == "ESE" ~ 112,
dPM == "SE" ~ 135,
dPM == "SSE" ~ 158,
dPM == "S" ~ 180,
dPM == "SSW" ~ 202,
dPM == "SW" ~ 225,
dPM == "WSW" ~ 248,
dPM == "W" ~ 270,
dPM == "WNW" ~ 292,
dPM == "NW" ~ 315,
dPM == "NNW" ~ 338),
specHumAM = (rhAM/100)*frame:::QSat(tAM, MSLP),
specHumPM = (rhPM/100)*frame:::QSat(tPM, MSLP))
# Create sequence
Temp <- clim[,c('tAM','tPM','tMin','tMax','Day')] %>% pivot_longer(c(tMin, tAM, tMax, tPM), names_to = "Time", values_to = "Temp") %>%
mutate(Time = case_when(Time == "tMin" ~ 3,
Time == "tAM" ~ 9,
Time == "tMax" ~ 12,
Time == "tPM" ~ 15))
Temp$Hour <- (Temp$Day * 24 - 24) + Temp$Time
RH <- clim[,c('rhAM','rhPM','Day')] %>% pivot_longer(c(rhAM, rhPM), names_to = "Time", values_to = "RH") %>%
mutate(Time = case_when(Time == "rhAM" ~ 9,
Time == "rhPM" ~ 15))
specHum <- clim[,c('specHumAM','specHumPM','Day')] %>% pivot_longer(c(specHumAM, specHumPM), names_to = "Time", values_to = "specHum") %>%
mutate(Time = case_when(Time == "specHumAM" ~ 9,
Time == "specHumPM" ~ 15))
Wind <- clim[,c('wAM','wPM','Day')] %>% pivot_longer(c(wAM, wPM), names_to = "Time", values_to = "Wind") %>%
mutate(Time = case_when(Time == "wAM" ~ 9,
Time == "wPM" ~ 15))
cloud <- clim[,c('cAM','cPM','Day')] %>% pivot_longer(c(cAM, cPM), names_to = "Time", values_to = "Cloud") %>%
mutate(Time = case_when(Time == "cAM" ~ 9,
Time == "cPM" ~ 15))
Direction <- clim[,c('dAM','dPM','Day')] %>% pivot_longer(c(dAM, dPM), names_to = "Time", values_to = "Direction") %>%
mutate(Time = case_when(Time == "dAM" ~ 9,
Time == "dPM" ~ 15))
Rain <- clim[,c('Rain', 'MSLP', 'Day')] %>%
mutate(Time = 9)
out <- left_join(Temp, cloud, by = c("Day", "Time"))
out <- left_join(out, RH, by = c("Day", "Time"))
out <- left_join(out, specHum, by = c("Day", "Time"))
out <- left_join(out, Wind, by = c("Day", "Time"))
out <- left_join(out, Rain, by = c("Day", "Time"))
out <- left_join(out, Direction, by = c("Day", "Time"))
out <- out %>%
select(Day, Hour, Rain, MSLP, Cloud, Temp, RH, specHum, Wind, Direction)
out <- out[order(out$Hour),]
return(out)
}
#' Formats hourly weather data for frame, with DFMC
#'
#' Function uses the 'Single differential equation model' of Matthews et al (2010)
#'
#' @param clim A dataset with the fields:
#' tAM (9am temp, degC)
#' tPM (3pm temp, degC)
#' tMin (Minimum daily temp, degC)
#' tMax (Maximum daily temp, degC)
#' rhAM (9am Relative humidity, %)
#' rhPM (3pm Relative humidity, %)
#' wAM (9am wind, km/h)
#' wPM (3pm wind, km/h)
#' cAM (Morning cloud, oktas)
#' cPM (Afternoon cloud, oktas)
#' MSLP (Mean Sea Level Pressure, hPa)
#' Rain (Daily rainfall, mm)
#' @param m The starting moisture (proportion ODW)
#' @param nelsonA Constant from Nelson FMC model
#' @param nelsonB Constant from Nelson FMC model
#' @param sigma Surface area: volume ratio of litter particles (m2/m3)
#' @param conLitter Litter conductance (m/s). Using value from Matthews spreadsheet rather than paper
#' @param dt Seconds per time step
#' @param LAI Leaf area index of the vegetation
#' @param WRF Wind reduction factor
#' @param rholitter Density of litter particles (kg/m3)
#' @param litterW Weight of litter (t/ha)
#' @param lat Latitude (degrees)
#' @param altitude (m)
#' @param slope (degrees)
#' @param rangeDir Cardinal direction of the ridgelines - either Nth/Sth (0) or west/east (270)
#' @param EPS Sensitivity of iteration
#' @param hCan Height of the tree canopy (m)
#' @param cardinal Simplify compass directions to cardinal (TRUE/FALSE)
#' @param slopeSD Standard deviation of the slope
#'
#' @export
frameWeather <- function(clim, m = 0.15, LAI = 3, WRF = 3, hCan = 20, rholitter = 550, litterW = 10,
lat = -31.89, altitude = 8, slope = 6.7, slopeSD = 4.9, rangeDir = 270, dt = 3600,
cardinal = FALSE, nelsonA = 5.2, nelsonB = -19, conLitter = 0.0006, sigma = 3000, EPS = 0.01) {
climDay <- tidyWeather(clim)
Rain <- climDay[,c('Hour', 'Rain')]
Direction <- climDay[,c('Hour', 'Direction')]
if (slope != 0) {
cat("Slope set to zero until issues resolved")
}
# Interpolate
Hour <- seq(from = 1, to = max(climDay$Hour), by = 1)
TempI <- splines::interpSpline( climDay$Hour, climDay$Temp )
MSLPI <- splines::interpSpline( climDay$Hour, climDay$MSLP, na.action = na.omit )
CloudI <- splines::interpSpline( climDay$Hour, climDay$Cloud, na.action = na.omit )
sRHI <- splines::interpSpline( climDay$Hour, climDay$specHum, na.action = na.omit )
WindI <- splines::interpSpline( climDay$Hour, climDay$Wind, na.action = na.omit )
Temp <- predict( TempI, Hour )$y
MSLP <- predict( MSLPI, Hour )$y
Cloud <- predict( CloudI, Hour )$y
sRH <- predict( sRHI, Hour )$y
Wind <- predict( WindI, Hour )$y
# Combine data
hemisphere <- if (lat<0) {
-1
} else {1}
out <- as.data.frame(Hour)
out$TempA <- Temp
# Adjust temp for canopy height. Model taken from Fig. 2 of
# Jucker, T., Hardwick, S.R., Both, S., Elias, D.M.O., Ewers, R.M., Milodowski, D.T., et al. (2018).
# Canopy structure and topography jointly constrain the microclimate of human-modified tropical landscapes. Glob. Chang. Biol., 24, 5243–5258.
out$Temp = (-(0.0009*exp(0.1195*(out$TempA-273.15)))*log(hCan)+1)*(out$TempA-273.15)+273.15
out$sRH <- sRH
out$MSLP <- MSLP
out$RHA <- (sRH / frame:::QSat(out$TempA, out$MSLP))*100
out$RH <- (sRH / frame:::QSat(out$Temp, out$MSLP))*100
out$Wind <- Wind / WRF
out$Cloud <- Cloud
# Add rain
out <- out %>%
left_join(Rain, by = "Hour")
out$Rain[is.na(out$Rain)] <- 0
# Add direction
# Optional set to cardinal directions
if (cardinal == TRUE) {
Direction <- Direction %>%
mutate(Direction = case_when(Direction == 22 ~ 0,
Direction == 45 ~ 0,
Direction == 68 ~ 90,
Direction == 112 ~ 90,
Direction == 135 ~ 90,
Direction == 158 ~ 180,
Direction == 202 ~ 180,
Direction == 225 ~ 180,
Direction == 248 ~ 270,
Direction == 292 ~ 270,
Direction == 315 ~ 270,
Direction == 338 ~ 0,
TRUE ~ Direction))
}
# First value
r <- 1
nr <- NA
while (is.na(nr)) {
nr <- Direction$Direction[r]
r <- r+1
}
out <- left_join(out, Direction, by = "Hour")
out$Direction[1] <- nr
for (r in 2:nrow(out)) {
if (is.na(out$Direction[r])) {
out$Direction[r] <- out$Direction[r-1]
}
}
# Calculate inputs & moisture
out <- out %>%
mutate(RH = pmax(0, pmin(100, RH)),
Wind = pmax(0, Wind),
Cloud = pmax(0, pmin(1, Cloud)),
LAI = LAI,
WRF = WRF,
Declination = hemisphere*23.45*cos((360*((Hour/24)-172)/365)*pi/180),
hourAngle = abs(12-((Hour/24)-floor(Hour/24))*24)*15,
solarAltitude = asin(cos(lat*pi/180)*cos(Declination*pi/180)*cos(hourAngle*pi/180)+sin(lat*pi/180)*sin(Declination*pi/180)),
#TerrainA = (cos((hourAngle-(rangeDir+90))*pi/180)*slope)*pi/180,
#TerrainB = (cos((hourAngle-(rangeDir-90))*pi/180)*slope)*pi/180,
#TerrainC = (cos((hourAngle-(rangeDir-90))*pi/180)*(slope+slopeSD))*pi/180,
#TerrainD = (cos((hourAngle-(rangeDir-90))*pi/180)*(slope+2*slopeSD))*pi/180,
# Set all terrain effects to 0 (flat) because model is drawing slopeSD from somewhere unknown
TerrainA = 0,
TerrainB = 0,
TerrainC = 0,
TerrainD = 0,
ZenithA = acos(sin(lat*pi/180)*sin(Declination*pi/180)+cos(lat*pi/180)*cos(Declination*pi/180)*cos(hourAngle*pi/180))-TerrainA,
ZenithB = acos(sin(lat*pi/180)*sin(Declination*pi/180)+cos(lat*pi/180)*cos(Declination*pi/180)*cos(hourAngle*pi/180))-TerrainB,
ZenithC = acos(sin(lat*pi/180)*sin(Declination*pi/180)+cos(lat*pi/180)*cos(Declination*pi/180)*cos(hourAngle*pi/180))-TerrainC,
ZenithD = acos(sin(lat*pi/180)*sin(Declination*pi/180)+cos(lat*pi/180)*cos(Declination*pi/180)*cos(hourAngle*pi/180))-TerrainD,
ShadeA = 1-exp(-(2/(pi*tan(pmax(0.01,1.570796-ZenithA))))*LAI),
ShadeB = 1-exp(-(2/(pi*tan(pmax(0.01,1.570796-ZenithB))))*LAI),
ShadeC = 1-exp(-(2/(pi*tan(pmax(0.01,1.570796-ZenithC))))*LAI),
ShadeD = 1-exp(-(2/(pi*tan(pmax(0.01,1.570796-ZenithD))))*LAI),
InsolationA = pmax(0,((1-ShadeA-Cloud)*1000*cos(ZenithA))),
InsolationB = pmax(0,((1-ShadeB-Cloud)*1000*cos(ZenithB))),
InsolationC = pmax(0,((1-ShadeC-Cloud)*1000*cos(ZenithC))),
InsolationD = pmax(0,((1-ShadeD-Cloud)*1000*cos(ZenithD))),
RainAdj = pmax(0,Rain-(0.001*(Rain/2)+0.08*LAI)),
Wetting =((10*pmin(RainAdj,(0.01*RainAdj+0.36*(litterW/5))))/litterW),
# Leaf temperatures at the soil, uses function HeatingModel in Matthews spreadsheet, Eq. 3 in Matthews et al (2010)
SoilA = Temp + (8.15 - 2.25 * exp(-0.6 * Wind) - 0.0312 * Temp +
(0.021 + (-0.04 + 0.0006 * Temp - 0.00000125 * Temp ^ 2) * exp(-0.6 * Wind)) * InsolationA) - 273.15, #Equation from spreadsheet
SoilB = Temp + (8.15 - 2.25 * exp(-0.6 * Wind) - 0.0312 * Temp +
(0.021 + (-0.04 + 0.0006 * Temp - 0.00000125 * Temp ^ 2) * exp(-0.6 * Wind)) * InsolationB) - 273.15,
SoilC = Temp + (8.15 - 2.25 * exp(-0.6 * Wind) - 0.0312 * Temp +
(0.021 + (-0.04 + 0.0006 * Temp - 0.00000125 * Temp ^ 2) * exp(-0.6 * Wind)) * InsolationC) - 273.15,
SoilD = Temp + (8.15 - 2.25 * exp(-0.6 * Wind) - 0.0312 * Temp +
(0.021 + (-0.04 + 0.0006 * Temp - 0.00000125 * Temp ^ 2) * exp(-0.6 * Wind)) * InsolationD) - 273.15)
for (t in out$Hour) {
if (t == 1) {
mA <- m
mB <- m
mC <- m
mD <- m
} else {
mA <- out$moistureA[t-1]
mB <- out$moistureB[t-1]
mC <- out$moistureC[t-1]
mD <- out$moistureD[t-1]
}
out$moistureA[t] <- pmax(0.01,(frame::simplefmc(m = mA,
tAir = out$Temp[t],
vAir = out$Wind[t],
pAir = out$MSLP[t],
rhAir = out$RH[t],
insolation = out$InsolationA[t],
rhoLitter = rholitter,
dt = dt,
nelsonA = nelsonA,
nelsonB = nelsonB,
conLitter = conLitter,
sigma = sigma,
EPS = EPS))) - max(0,0.5*(mA-1)) + out$Wetting[t]
out$moistureB[t] <- pmax(0.01,(frame::simplefmc(m = mB,
tAir = out$Temp[t],
vAir = out$Wind[t],
pAir = out$MSLP[t],
rhAir = out$RH[t],
insolation = out$InsolationB[t],
rhoLitter = rholitter,
dt = dt,
nelsonA = nelsonA,
nelsonB = nelsonB,
conLitter = conLitter,
sigma = sigma,
EPS = EPS))) - max(0,0.5*(mB-1)) + out$Wetting[t]
out$moistureC[t] <- pmax(0.01,(frame::simplefmc(m = mC,
tAir = out$Temp[t],
vAir = out$Wind[t],
pAir = out$MSLP[t],
rhAir = out$RH[t],
insolation = out$InsolationC[t],
rhoLitter = rholitter,
dt = dt,
nelsonA = nelsonA,
nelsonB = nelsonB,
conLitter = conLitter,
sigma = sigma,
EPS = EPS))) - max(0,0.5*(mC-1)) + out$Wetting[t]
out$moistureD[t] <- pmax(0.01,(frame::simplefmc(m = mD,
tAir = out$Temp[t],
vAir = out$Wind[t],
pAir = out$MSLP[t],
rhAir = out$RH[t],
insolation = out$InsolationD[t],
rhoLitter = rholitter,
dt = dt,
nelsonA = nelsonA,
nelsonB = nelsonB,
conLitter = conLitter,
sigma = sigma,
EPS = EPS))) - max(0,0.5*(mD-1)) + out$Wetting[t]
}
out <- out %>%
mutate(Temp = Temp - 273.15,
TempA = TempA - 273.15,
Wind = Wind * 3.6,
MSLP = MSLP/100,
wetBulb = TempA*atan(0.151977*(RHA+8.313659)^0.5)+atan(TempA+RHA)-atan(RHA-1.676331)+0.00391838*RHA^(3/2)*atan(0.023101*RHA)-4.686035,
cgStrikes = 2.708*10^-46*exp(3.863*wetBulb))
return(out)
}
#' Internal function for climDynamics
#'
#' @param a
#'
#' @return dataframe
#' @export
#'
parClim <- function(a) {
FloraA <- filter(Flora, record == a)
StructureA <- filter(Structure, record == a)
base.params <- suppressWarnings(frame::buildParams(StructureA, FloraA, default.species.params, a,
fLine = 1, slope = 0, temp = 30, dfmc = 0.05, wind = 10))
hCan <- max(as.numeric(FloraA$top), na.rm = TRUE)
LAI <- LAIcomm(base.params, yu = hCan, yl = 0) * as.numeric(lAngleAccounting) * ((as.numeric(!is.na(StructureA$NS))+as.numeric(!is.na(StructureA$E))+as.numeric(!is.na(StructureA$M))+as.numeric(!is.na(StructureA$C)))/4)
WRF <- windReduction(base.params, test = 1.2)
litterW <- as.numeric(max(FloraA$weight, na.rm = TRUE))
out <- frame::frameWeather(clim = clim, m, LAI, WRF, hCan, rholitter, litterW,
lat, slope, slopeSD, rangeDir, cardinal) %>%
mutate(Record = a,
site = StructureA$site[1])
return(out)
}
#' Models input weather parameters from a climate dataset,
#' each age is modelled on a separate core
#'
#' @param fireDat
#' @param clim
#' @param m
#' @param rholitter
#' @param lat Latitude (degrees)
#' @param slope (degrees)
#' @param slopeSD Standard deviation of the slope
#' @param rangeDir Cardinal direction of the ridgelines - either Nth/Sth (0) or west/east (270)
#' @param cardinal
#' @param freeCores Number of cores to leave unused by the model
#' @param lAngleAccounting Set to TRUE to automatically modify leaf angle by number of strata
#'
#' @return dataframe
#' @export
#'
climDynamics <- function(fireDat, clim, m = 0.15, slope = 0, slopeSD = 0, rholitter = 550,
lat = -34.95, lAngleAccounting = TRUE, rangeDir = 270, cardinal = TRUE, freeCores = 1){
# 1. Compile inputs
Flora <- fireDat[[1]]
Structure <- fireDat[[2]]
default.species.params <- fireDat[[3]]
r <- unique(Flora$record)
Out <- data.frame()
# 2. Create a cluster of cores with replicated R on each
nCores <- max(parallel::detectCores() - freeCores,1)
cl <- parallel::makeCluster(nCores)
# 3. Load the packages
parallel::clusterEvalQ(cl,
{ library(dplyr)
library(tidyr)
library(frame)
library(assertthat)
library(extraDistr)})
# 4. Load the inputs
parallel::clusterExport(cl,varlist=c('Flora', 'Structure', 'default.species.params', 'clim', 'm',
'slope', 'slopeSD', 'rholitter', 'lat', 'lAngleAccounting', 'rangeDir', 'cardinal'), environment())
# 5. Send each rep to a different core to be processed
system.time(out <- parallel::parLapply(cl, r, parClim))
parallel::stopCluster(cl)
for (n in 1:length(r)) {
Na <- as.data.frame(out[[n]])
Out <- rbind(Out, Na)
}
return(Out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.