R/NicheMapRcode.R
In ilyamaclean/microclimc: Below or above canopy microclimate modelling with R
Defines functions
runwithNMR
runmodelS
.runmodelsnow
.snowtempf
tleafS
.windprofile
.leafabs2
runNMR
Documented in
.leafabs2
runmodelS
.runmodelsnow
runNMR
runwithNMR
.snowtempf
tleafS
.windprofile
#' Runs NicheMapR model
#'
#' @description Wrapper function for running NicheMapR to derive soil moistures
#' @param climdata a data.frame of hourly weather variables in same format as [microclimc::weather()].
#' Column obs_time must give times in GMT/UTC.
#' @param prec a vector of daily or hourly precipitation values (mm).
#' @param lat latitude (decimal degrees, positive in northern hemisphere).
#' @param long longitude (decimal degrees, negative west of Greenwich meridian).
#' @param Usrhyt Local height (m) at which air temperature, wind speed and humidity are to be computed (see details).
#' @param Veghyt At of vegetation canopy (see details).
#' @param Refhyt Reference height (m) at which air input climate variables are measured.
#' @param PAI single value or vector of hourly or daily values of total plant area per unit ground area of canopy.
#' @param LOR Campbell leaf angle distrubution coefficient
#' @param pLAI fraction of PAI that is green vegetation.
#' @param clump clumpiness factor for canopy (0 to 1, 0 = even)
#' @param REFL A single numeric value of ground reflectivity of shortwave radiation.
#' @param LREFL A single numeric value of average canopy reflectivity of shortwave radiation.
#' @param SLE Thermal emissivity of ground
#' @param DEP Soil depths at which calculations are to be made (cm), must be 10 values
#' starting from 0, and more closely spaced near the surface.
#' @param ALTT elevation (m) of location.
#' @param SLOPE slope of location (decimal degrees).
#' @param ASPECT aspect of location (decimal degrees, 0 = north).
#' @param ERR Integrator error tolerance for soil temperature calculations.
#' @param soiltype soil type as for [microclimc::soilparams()]. Set to NA to use soil parameters
#' specified below.
#' @param PE Air entry potentials (J/kg) (19 values descending through soil for specified soil
#' nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param KS Saturated conductivities (kg s/m^3) (19 values descending through soil for
#' specified soil nodes in parameter DEP and points half way between). Ignored if soil type
#' is given.
#' @param BB Campbell's soil 'b' parameter (19 values descending through soil for specified
#' soil nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param BD Soil bulk density (Mg/m^3) (19 values descending through soil for specified
#' soil nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param DD Soil density (Mg/m^3) (19 values descending through soil for specified
#' soil nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param cap Is organic cap present on soil surface? (1 = Yes, 0 = No).
#' @param hori Horizon angles (degrees), from 0 degrees azimuth (north) clockwise in 10 degree intervals.
#' @param maxpool Max depth for water pooling on the surface (mm), to account for runoff.
#' @param rainmult Rain multiplier for surface soil moisture (used to induce runon).
#' @param SoilMoist_Init Initial volumetric soil water content at each soil node (m^3/m^3)
#' @return a list with the following components:
#' @return `metout` a data.frame of microclimate variables for each hour at height `Usrhyt`.
#' @return `soiltemps` a data.frame of hourly soil temperatures for each depth node.
#' @return `soilmoist` a data.frame of hourly soil volumetric water content for each depth node.
#' @return `snowtemp` if snow present, a data,frame of snow temperature (°C), at each of the potential 8 snow layers (see details).
#' 0 if snow not present.
#' @return `plant` a data.frame of plant transpiration rates (g/m^2/hr), leaf water potentials
#' (J/kg) and root water potential (J/kg) at each of the 10 specified depths.
#' @return `nmrout` full NicheMapR output
#' @details Requires NicheMapR: devtools::install_github('mrke/NicheMapR'). NicheMapR is an integrated
#' soil moisture and temperature model that treats the vegetation as a single layer (https://mrke.github.io/).
#' This wrapper function, runs the NicheMapR::microclimate function with reduced parameter
#' inputs, by specifying sensible default values for other parameters where it is unlikely that these
#' would be known or explicitely measured at the site. The degree of canopy shading is worked out
#' explicitely from `PAI` at each hourly time interval by estimating canopy tranmission of direct
#' and diffuse radiation, and by adjusting the amount of radiation that would be absorbed
#' at the ground surface for a given slope and aspect. Roughness lengths and zero plane displacement
#' heights are limited to <2m in NicheMapR, so if `Veghyt` > 2, it is set to 2m.
#' This is adequate for computing soil moisture and total evapotranspiration, but will give false
#' air temperatures for heights below canopy (see [runwithNMR()]. If `soiltype` is given, the subsequent
#' soil parameters are computed from soil type. If `soiltype` is NA, soil paramaters must be specified.
#'
#' @author Ilya Maclean (i.m.d.maclean@exeter.ac.uk) and Urtzi Urzelai (urtzi.enriquez@gmail.com)
#' @import microctools
#' @import NicheMapR
#' @export
#' @examples
#' # Run NicheMapR with default parameters and inbuilt weather datasets
#' library(NicheMapR)
#' # Plot air temperatures
#' microout<-runNMR(weather,dailyprecip,50.2178,-5.32656,0.05,0.1,PAI=1)
#' metout <- microout$metout
#' tmn <- min(metout$TALOC,metout$TAREF)
#' tmx <- max(metout$TALOC,metout$TAREF)
#' dday <- metout$DOY+metout$TIME/1440 # Decimal day
#' plot(metout$TALOC~dday, type="l", col = "red", ylim=c(tmn,tmx), ylab = "Temperature")
#' par(new = T)
#' plot(metout$TAREF~dday, type="l", col = "blue", ylim=c(tmn,tmx), ylab = "", xlab = "")
#' # Plot soil temperatures
#' soiltemp<-microout$soiltemps
#' tmn <- min(soiltemp$D0cm,soiltemp$D30cm)
#' tmx <- max(soiltemp$D0cm,soiltemp$D30cm)
#' plot(soiltemp$D0cm~dday, type="l", col = "red", ylim=c(tmn,tmx), ylab = "Temperature")
#' par(new = T)
#' plot(soiltemp$D30cm~dday, type="l", col = "blue", ylim=c(tmn,tmx), ylab = "", xlab = "")
#' # Plot soil moistures
#' soilm <- microout$soilmoist
#' mmn <- min(soilm$WC2.5cm,soilm$WC30cm)
#' mmx <- max(soilm$WC2.5cm,soilm$WC30cm)
#' plot(soilm$WC2.5cm~dday, type="l", col = "red", ylim=c(mmn,mmx), ylab = "Soil moisture")
#' par(new = T)
#' plot(soilm$WC30cm~dday, type="l", col = "blue", ylim=c(mmn,mmx), ylab = "", xlab = "")
runNMR <- function(climdata, prec, lat, long, Usrhyt, Veghyt, Refhyt = 2, PAI = 3, LOR = 1,
pLAI = 0.8, clump = 0, REFL = 0.15, LREFL = 0.4, SLE = 0.95,
DEP = c(0,2.5,5,10,15,20,30,50,100,200), ALTT = 0, SLOPE = 0, ASPECT = 0,
ERR = 1.5, soiltype = "Loam", PE = NA, KS = NA, BB = NA, BD = NA, DD = NA,
cap = 1, hori = rep(0,36), maxpool = 1000, rainmult = 1,
SoilMoist_Init = c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3),
animal = FALSE) {
# Internal functions
.zeroplanedis<-function(h,pai) {
pai[pai<0.001]<-0.001
d<-(1-(1-exp(-sqrt(7.5*pai)))/sqrt(7.5*pai))*h
d
}
#' Calculate roughness length
.roughlength<-function(h,pai,d) {
Be<-sqrt(0.003+(0.2*pai)/2)
zm<-(h-d)*exp(-0.4/Be)
zm
}
# Hourly to daily
.dmaxmin <- function(x, fun) {
dx <- t(matrix(x, nrow = 24))
apply(dx, 1, fun)
}
# Get soil parameters
.getsoilparams<-function(soiltype) {
# Get soil parameters based on soil type
sel<-which(soilparams$Soil.type==soiltype)
BulkDensity<-soilparams$rho[sel]
SatWater <- soilparams$Smax[sel] # volumetric water content at saturation (0.1 bar matric potential) (m3/m3)
Clay <- CampNormTbl9_1$Clay[sel]*100 # clay con
KS<-rep(CampNormTbl9_1$Ks[sel],19)
BB<-rep(CampNormTbl9_1$b[sel],19)
PE<-rep(CampNormTbl9_1$airentry[sel],19)
return(list(BulkDensity=BulkDensity,SatWater=SatWater,Clay=Clay,KS=KS,BB=BB,PE=PE))
}
# Time and location parameters
ndays<-length(climdata$temp)/24
tme<-as.POSIXlt(climdata$obs_time,tz="UTC")
doy<-.dmaxmin(tme$yday+1,mean)
idayst <- 1 # start day (legacy parameter)
ida <- ndays # end day (legacy parameter)
HEMIS <- ifelse(lat < 0, 2, 1)
ALAT<-abs(trunc(lat))
AMINUT<-(abs(lat)-ALAT)*60
ALONG<- abs(trunc(long))
ALMINT<-(abs(long)-ALONG)*60
ALREF<-0
EC <- 0.0167238
# Air and wind vertical profile parameters
D0<-.zeroplanedis(Veghyt,mean(PAI))
RUF<-.roughlength(Veghyt,mean(PAI),D0)
D0<-0
ZH<-0.2*RUF
# Terrain and shading parameters
VIEWF<-1-sum(sin(hori*pi/180))/length(hori)
# Soil properties
Numtyps <- 2 # number of soil types
Nodes <- matrix(data = 0, nrow = 10, ncol = ndays) # array of all possible soil nodes
Nodes[1,1:ndays] <- 3 # deepest node for first substrate type
Nodes[2,1:ndays] <- 9 # deepest node for second substrate type
# Time varying environmental data
tannul <- mean(climdata$temp)
# Shade parameter
grasshade<-ifelse(Veghyt<0.5,1,0) # VT
if (length(PAI) == 1) {
PAI<-rep(PAI,ndays)
} else if (length(PAI) == ndays*24) {
PAI<-matrix(PAI,ncol=24,byrow=T)
PAI<-apply(PAI,1,mean)
}
if (length(PAI) != ndays) stop("PAI must be a single value or hourly/daily \n")
PAIc<-PAI^(1/(1-clump))
MINSHADES<-((1-exp(-PAIc))+clump)*100
MAXSHADES<-rep(100,ndays)
# Modal times of air temp, wind, humidity and cloud cover
TIMINS <- c(0, 0, 1, 1)
TIMAXS <- c(1, 1, 0, 0)
# Spinup
soilinit <- rep(tannul, 20)
spinup <- 1
# Decide whether to run snow model
if (length(prec) == ndays) {
RAINhr<-rep(prec/24,each=24)
} else RAINhr<-prec
snowmodel<-max(ifelse(climdata$temp<0,1,0)*ifelse(RAINhr>0,0,1))
densfun <- c(0.5979, 0.2178, 0.001, 0.0038)
intercept <- mean(MINSHADES) / 100 * 0.3
# Decide whether to run snowmodel
microinput<-c(ndays, RUF, ERR, Usrhyt, Refhyt, Numtyps, Z01 = 0, Z02 = 0, ZH1 = 0,
ZH2 = 0, idayst = 1, ida, HEMIS, ALAT, AMINUT, ALONG, ALMINT, ALREF,
SLOPE, ASPECT, ALTT, CMH2O = 1, microdaily = 1, tannul, EC, VIEWF,
snowtemp = 1.5, snowdens = 0.375, snowmelt = 0.9, undercatch = 1,
rainmult, runshade = 1, runmoist = 1, maxpool, evenrain = 1, snowmodel,
rainmelt = 0.0125, writecsv = 0, densfun, hourly = 1, rainhourly = 0,
lamb = 0, IUV = 0, RW = 2.5e+10, PC = -1500, RL = 2e+6, SP = 10, R1 = 0.001,
IM = 1e-06, MAXCOUNT = 500, IR = 0, message = 0, fail = 8760, snowcond = 0,
intercept, grasshade, solonly = 0, ZH, D0, TIMAXS, TIMINS, spinup,0,360)
tides <- matrix(data = 0, nrow = 24*ndays, ncol = 3)
SLES <- rep(SLE, ndays)
# Weather variables
TAIRhr<-climdata$temp
RHhr<-climdata$relhum
WNhr<-climdata$windspeed
PRESShr<-climdata$pres*1000
ea<-satvap(TAIRhr)*(RHhr/100)
eo<-1.24*(10*ea/(TAIRhr+273.15))^(1/7)
CLDhr<-((climdata$skyem-eo)/(1-eo))*100
CLDhr[CLDhr<0]<-0
CLDhr[CLDhr>100]<-100
TMAXX<-.dmaxmin(TAIRhr,max)
TMINN<-.dmaxmin(TAIRhr,min)
CCMAXX<-.dmaxmin(CLDhr,max)
CCMINN<-.dmaxmin(CLDhr,min)
RHMAXX<-.dmaxmin(RHhr,max)
RHMINN<-.dmaxmin(RHhr,min)
WNMAXX<-.dmaxmin(WNhr,max)
WNMINN<-.dmaxmin(WNhr,min)
PRESS<-.dmaxmin(PRESShr,min)
SOLRhr<-climdata$swrad*VIEWF
sb<-5.67*10^-8
IRDhr<-climdata$skyem*sb*(climdata$temp+273.15)^4
lt<-tme$hour+tme$min/60+tme$sec/3600
jd<-jday(tme=tme)
sa<-solalt(lt,lat,long,jd,0)
ZENhr<-90-sa
ZENhr[ZENhr>90]<-90
REFLS <- rep(REFL, ndays) # set up vector of soil reflectances for each day
# Soil surface wetness
RAIND<-.dmaxmin(RAINhr,max)
PCTWET<-ifelse(RAIND>0,90,0)
# Aerosol extinction coefficient profile
optdep.summer<-as.data.frame(rungads(lat,long,1,0))
optdep.winter<-as.data.frame(rungads(lat,long,1,0))
optdep<-cbind(optdep.winter[,1],rowMeans(cbind(optdep.summer[,2],optdep.winter[,2])))
optdep<-as.data.frame(optdep)
colnames(optdep)<-c("LAMBDA","OPTDEPTH")
a<-lm(OPTDEPTH~poly(LAMBDA,6,raw=TRUE),data=optdep)
LAMBDA<-c(290,295,300,305,310,315,320,330,340,350,360,370,380,390,400,420,440,460,480,
500,520,540,560,580,600,620,640,660,680,700,720,740,760,780,800,820,840,860,
880,900,920,940,960,980,1000,1020,1080,1100,1120,1140,1160,1180,1200,1220,
1240,1260,1280,1300,1320,1380,1400,1420,1440,1460,1480,1500,1540,1580,1600,
1620,1640,1660,1700,1720,1780,1800,1860,1900,1950,2000,2020,2050,2100,2120,
2150,2200,2260,2300,2320,2350,2380,2400,2420,2450,2490,2500,2600,2700,2800,
2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000)
TAI<-predict(a,data.frame(LAMBDA))
# Soil properties (mineral component only)
Thcond <- 2.5 # soil minerals thermal conductivity (W/mC)
SpecHeat <- 870 # soil minerals specific heat (J/kg-K)
if (class(DD) == "logical") {
Density<-2.560 # soil minerals density (Mg/m3)
} else Density<-mean(DD)
SP<-.getsoilparams(soiltype)
if (class(BD) == "logical") {
BulkDensity <- SP$BulkDensity # soil bulk density (kg/m3)
} else BulkDensity<-mean(BD)
SatWater<-SP$SatWater # volumetric water content at saturation (0.1 bar matric potential) (m3/m3)
Clay <- SP$Clay # clay content for matric potential calculations (%)
# Create matrix of soil properties
soilprops <- matrix(data = 0, nrow = 10, ncol = 6) # create an empty soil properties matrix
soilprops[1:2,1] <- BulkDensity # insert soil bulk density to profile 1
soilprops[1:2,2] <- SatWater # insert saturated water content to profile 1
soilprops[1:2,3] <- Thcond # insert thermal conductivity to profile 1
soilprops[1:2,4] <- SpecHeat # insert specific heat to profile 1
soilprops[1:2,5] <- Density # insert mineral density to profile 1
if(cap == 1) { # insert thermal conductivity to profile 1, and see if 'organic cap' added on top
soilprops[1,3]<-0.2 # mineral thermal conductivity
soilprops[1,4]<-1920
}
soilinit <- rep(tannul,20) # make inital soil temps equal to mean annual
if (class(PE) == "logical") PE<-SP$PE
if (class(KS) == "logical") KS<-SP$KS
if (class(BB) == "logical") BB<-SP$BB
if (class(BD) == "logical") BD<-rep(BulkDensity,19)
if (class(DD) == "logical") DD<-rep(Density,19)
# Initial soil moistures
moists <- matrix(nrow=10, ncol = ndays, data = 0) # set up an empty vector for soil moisture values through time
moists[1:10,] <- SoilMoist_Init # insert inital soil moisture
# Calculate monthly mean rainfall
RAINFALL<-.dmaxmin(RAINhr, sum)
tannulrun <- rep(tannul,ndays) # deep soil temperature (2m) (deg C)
# Root density at each node:
L<-c(0,0,8.2,8,7.8,7.4,7.1,6.4,5.8,4.8,4,1.8,0.9,0.6,0.8,0.4,0.4,0,0)*10000
LAI<-pLAI*PAI
micro<-list(microinput = microinput, tides=tides, doy = doy, SLES = SLES, DEP = DEP,
Nodes=Nodes, MAXSHADES = MAXSHADES, MINSHADES = MINSHADES, TMAXX = TMAXX,
TMINN = TMINN, RHMAXX = RHMAXX, RHMINN = RHMINN, CCMAXX = CCMAXX, CCMINN = CCMINN,
WNMAXX = WNMAXX, WNMINN = WNMINN, TAIRhr = TAIRhr, RHhr = RHhr, WNhr = WNhr,
CLDhr = CLDhr, SOLRhr = SOLRhr, RAINhr = RAINhr, ZENhr = ZENhr, IRDhr = IRDhr,
REFLS = REFLS, PCTWET = PCTWET, soilinit = soilinit, hori = hori,
TAI = TAI, soilprops = soilprops, moists = moists, RAINFALL = RAINFALL,
tannulrun = tannulrun, PE = PE, KS = KS, BB = BB, BD = BD, DD = DD,
L = L, LAI = LAI)
# Run model
microut<-microclimate(micro)
if(animal==FALSE){
metout<-as.data.frame(microut$metout)
soil<-as.data.frame(microut$soil)
soilmoist<-as.data.frame(microut$soilmoist)
plant<-as.data.frame(microut$plant)
if (snowmodel == 1) {
snow <- as.data.frame(microut$sunsnow)
} else snow <- 0
return(list(metout=metout,soiltemps=soil,soilmoist=soilmoist,snowtemp=snow,plant=plant,nmrout=microut))
}
if(animal == TRUE){
metout<-microut$metout
shadmet<-microut$shadmet
soil<-microut$soil
shadsoil<-microut$shadsoil
soilmoist<-microut$soilmoist
shadmoist <- microut$shadmoist
humid <- microut$humid
shadhumid <- microut$shadhumid
soilpot <- microut$soilpot
shadpot <- microut$shadpot
plant <- microut$plant
shadplant <- microut$shadplant
if (snowmodel == 1) {
snow <- microut$sunsnow
shdsnow <- microut$shdsnow
} else snow <- 0
drlam <- microut$drlam
drrlam <- microut$drrlam
srlam <- microut$srlam
timeinterval = 365
days <- rep(seq(1, timeinterval * nyears), 24)
days <- days[order(days)]
dates <- days + metout[, 2]/60/24 - 1
dates2 <- seq(1, timeinterval * nyears)
if (snowmodel == 1) {
nmrout_full <- list(soil = soil, shadsoil = shadsoil,
metout = metout, shadmet = shadmet, soilmoist = soilmoist,
shadmoist = shadmoist, humid = humid, shadhumid = shadhumid,
soilpot = soilpot, shadpot = shadpot, sunsnow = snow,
shdsnow = shdsnow, plant = plant, shadplant = shadplant,
RAINFALL = RAINFALL, ndays = ndays, elev = ALTT,
REFL = REFL[1], longlat = c(x[1], x[2]), nyears = nyears,
timeinterval = timeinterval, minshade = MINSHADES,
maxshade = MAXSHADES, DEP = DEP, drlam = drlam,
drrlam = drrlam, srlam = srlam, dates = dates,
dates2 = dates2)
}
else {
nmrout_full <- list(soil = soil, shadsoil = shadsoil,
metout = metout, shadmet = shadmet, soilmoist = soilmoist,
shadmoist = shadmoist, humid = humid, shadhumid = shadhumid,
soilpot = soilpot, shadpot = shadpot, plant = plant,
shadplant = shadplant, RAINFALL = RAINFALL,
ndays = ndays, elev = ALTT, REFL = REFL[1],
longlat = c(x[1], x[2]), nyears = nyears, timeinterval = timeinterval,
minshade = MINSHADES, maxshade = MAXSHADES,
DEP = DEP, drlam = drlam, drrlam = drrlam,
srlam = srlam, dates = dates, dates2 = dates2)
}
metout <- data.frame(metout)
soil <- data.frame(soil)
soilmoist <- data.frame(soilmoist)
snow <- data.frame(snow)
plant <- data.frame(plant)
return(list(metout=metout,soiltemps=soil,soilmoist=soilmoist,snowtemp=snow,plant=plant,nmrout=nmrout_full))
}
}
#' Internal function for calculating lead absorbed radiation on vector
.leafabs2 <-function(Rsw, tme, tair, tground, lat, long, PAIt, PAIu, pLAI, x, refls, refw, refg, vegem, skyem, dp,
merid = round(long/15, 0) * 15, dst = 0, clump = 0) {
jd<-jday(tme = tme)
lt<-tme$hour+tme$min/60+tme$sec/3600
sa<-solalt(lt,lat,long,jd,merid=merid,dst=dst)
sa2<-ifelse(sa<5,5,sa)
ref<-pLAI*refls+(1-pLAI)*refw
sunl<-psunlit(PAIu,x,sa,clump)
mul<-radmult(x,sa2)
aRsw <- (1-ref) * cansw(Rsw,dp,jd,lt,lat,long,PAIu,x,ref,merid=merid,dst=dst,clump=clump)
aGround <- refg * cansw(Rsw,dp,jd,lt,lat,long,PAIt,x,ref,merid=merid,dst=dst,clump=clump)
aGround <- (1-ref) * cansw(aGround,dp,jd,lt,lat,long,PAIt,x,ref,merid=merid,dst=dst,clump=clump)
aRsw <- sunl * mul * aRsw + (1-sunl) * aRsw + aGround
aRlw <- canlw(tair, PAIu, 1-vegem, skyem = skyem, clump = clump)$lwabs
return(aRsw+aRlw)
}
#' Internal function for calculating lead absorbed radiation on vector
.windprofile <- function(ui, zi, zo, a = 2, hgt, PAI, psi_m = 0, hgtg = 0.05 * hgt, zm0 = 0.004) {
d<-zeroplanedis(hgt,PAI)
zm<-roughlength(hgt,PAI,zm0)
ln2<-suppressWarnings(log((zi-d)/zm)+psi_m)
ln2[ln2<0.55]<-0.55
uf<-(0.4*ui)/ln2
if (zo >= hgt) {
ln2<-suppressWarnings(log((zo-d)/zm)+psi_m)
ln2[ln2<0.55]<-0.55
uo<-(uf/0.4)*ln2
} else {
ln2<-suppressWarnings(log((hgt-d)/zm)+psi_m)
ln2[ln2<0.55]<-0.55
uh<-(uf/0.4)*ln2
if (zi >= (0.1*hgt)) {
uo<-uh*exp(a*(zo/hgt)-1)
} else {
uo<-uh*exp(a*((0.1*hgt)/hgt)-1)
zmg<-0.1*hgtg
ln1<-log((0.1*hgt)/zmg)+psi_m
uf<-0.4*(uo/ln1)
ln2<-log(zo/zmg)+psi_m
uo<-(uf/0.4)*ln2
}
}
uo
}
#' Steady-state leaf and air temperature
#' @description Rapid method for calculating steady-state leaf and air temperature at one height
#' @param tair air temperature at reference height (deg C)
#' @param tground ground surface temperature (as returned by e.g. [runNMR()])
#' @param relhum relative humidity at reference height (percentage)
#' @param pk atmospheric pressure (kPa)
#' @param theta volumetric water content of upper most soil layer in current time step (m^3 / m^3)
#' as returned by e.g. [runNMR()]
#' @param gtt molar conductivity from leaf height to reference height as returned by [gturb()]
#' and [gcanopy()]
#' @param gt0 molar conductivity from leaf height to ground as returned by [gcanopy()]
#' @param gha boundary layer conductance of leaf as returned by [gforcedfree()]
#' @param gL combined boundary layer and leaf-air conductance given by 1/(1/gha+1/(uz*ph))
#' @param gv combined boundary layer and stomatal conductance of leaf
#' @param Rabs radiation absorbed by leaf as returned by e.g. [leafabs()]
#' @param soilb Shape parameter for Campbell soil model (dimensionless, > 1) as returned by
#' [soilinit()]
#' @param Psie Soil matric potential (J / m^3) as returned by [soilinit()]
#' @param Smax Volumetric water content at saturation (m3 / m3)
#' @param surfwet proportion of leaf surface acting as free water surface
#' @param leafdens Total one sided leaf area per m^3 at desired height
#' @return a list of the following:
#' @return `tleaf` leaf temperature
#' @return `tn` air temperature
#' @return `rh` relative humidity
#' @import microctools
#' @export
tleafS <- function(tair, tground, relhum, pk, theta, gtt, gt0, gha, gv, gL, Rabs, vegem, soilb,
Psie, Smax, surfwet, leafdens) {
###
cp<-cpair(tair)
# Air temperature expressed as leaf temperature
aL<-(gtt*(tair+273.15)+gt0*(tground+273.15))/(gtt+gt0)
bL<-(leafdens*gL)/(gtt+gt0)
# Vapour pressures
es<-satvap(tair)
eref <- (relhum/100)*es
rhsoil<-soilrh(theta,soilb,Psie,Smax,tground)
esoil<-rhsoil*satvap(tground)
# add a small correction to improve delta estimate
sb<-5.67*10^-8
Rnet<-Rabs-0.97*sb*(tair+273.15)^4
tle<-tair+(0.5*Rnet)/(cp*gha)
wgt1<-ifelse(leafdens<1,leafdens/2,1-0.5/leafdens)
wgt2<-1-wgt1
tapprox<-wgt1*tle+wgt2*tair
delta <- 4098*(0.6108*exp(17.27*tapprox/(tapprox+237.3)))/(tapprox+237.3)^2
ae<-(gtt*eref+gt0*esoil+gv*es)/(gtt+gt0+gv)
be<-(gv*delta)/(gtt+gt0+gv)
# Sensible heat
bH<-gha*cp
# Latent heat
lambda <- (-42.575*tair+44994)
aX<-((lambda*gv)/pk)*(surfwet*es-ae)
bX<-((lambda*gv)/pk)*(surfwet*delta-be)
aX[aX<0]<-0
bX[bX<0]<-0
# Emmited radiation
aR<-sb*vegem*aL^4
bR<-4*vegem*sb*(aL^3*bL+(tair+273.15)^3)
# Leaf temperature
dTL <- (Rabs-aR-aX)/(1+bR+bX+bH)
# tz pass 1
tn<-aL-273.15+bL*dTL
tleaf<-tn+dTL
# new vapour pressure
eanew<-ae+be*dTL
eanew[eanew<0.01]<-0.01
tmin<-dewpoint(eanew,tn,ice = TRUE)
esnew<-satvap(tn)
eanew<-ifelse(eanew>esnew,esnew,eanew)
rh<-(eanew/esnew)*100
# Set both tair and tleaf so as not to drop below dewpoint
tleaf<-ifelse(tleaf<tmin,tmin,tleaf)
tn<-ifelse(tn<tmin,tmin,tn)
tmax<-ifelse(tn+20<80,tn+20,80)
tleaf<-ifelse(tleaf>tmax,tmax,tleaf)
# cap upper limits of both tair and tleaf
tmx<-pmax(tground+5,tair+5)
tn<-ifelse(tn>tmx,tmx,tn)
tmx<-pmax(tground+20,tair+20)
tleaf<-ifelse(tleaf>tmx,tmx,tleaf)
# cap lower limits
tmn<-tair-7
tn<-ifelse(tn<tmn,tmn,tn)
tmn<-tair-20
tleaf<-ifelse(tleaf<tmn,tmn,tleaf)
return(list(tleaf=tleaf,tn=tn,rh=rh))
}
#' Internation function for computing snow temperature
.snowtempf <- function(snowdepth,snow,reqhgt) {
ndepths<-rev(c(0,2.5,5,10,20,50,100,200,300))/100
ha<-ndepths-reqhgt
selb<-which(ha<=0)[1]# layer below
sela<-selb-1 # layer below or equal
sm<-abs(reqhgt-ndepths[sela])+abs(reqhgt-ndepths[selb])
wgt1<-1-abs(reqhgt-ndepths[sela])/sm # weight above
wgt2<-1-abs(reqhgt-ndepths[selb])/sm # weight below
hs<-which(ndepths[sela]>(snowdepth/100))
wgt1<-rep(wgt1,length(snowdepth))
wgt2<-rep(wgt2,length(snowdepth))
wgt1[hs]<-0
wgt2[hs]<-1
stemp<-wgt1*snow[,sela+2]+wgt2*snow[,selb+2]
sel<-which(snowdepth<reqhgt)
stemp[sel]<-0
stemp
}
#' Internal function for running model with snow
.runmodelsnow <- function(climdata, vegp, soilp, nmrout, reqhgt, lat, long, metopen = TRUE, windhgt = 2) {
# Snow
snow<-nmrout$snow
if (class(snow)[1] == "data.frame") {
snowtemp<-snow$SN1
} else snowtemp<-rep(0,length(L))
metout<-nmrout$metout
snowdep<-metout$SNOWDEP/100
# (1) Unpack variables
tme<-as.POSIXlt(climdata$obs_time)
tair<-climdata$temp
relhum<-climdata$relhum
pk<-climdata$pres
u<-climdata$windspeed
hgt<-vegp$hgt
# Esimate PAIu
PAIt<-apply(vegp$PAI,2,sum)
PAIt[PAIt<0.001]<-0.001
LAI<-(vegp$pLAI*vegp$PAI)
LAI<-apply(LAI,2,mean)
pLAI<-LAI/PAIt
pLAI[is.na(pLAI)]<-0.8
m<-length(vegp$iw)
z<-c((1:m)-0.5)/m*hgt
# Esimate PAI above point and PAI of layer
if (reqhgt < hgt) {
sel<-which(z>reqhgt)
if (length(sel) > 1) {
wgt1<-abs(z[sel[1]]-reqhgt)
wgt2<-abs(z[sel[1]-1]-reqhgt)
sel2<-c(sel[1]-1,sel)
PAIu1<-vegp$PAI[sel,]
if (length(sel) > 1) PAIu1<-apply(PAIu1,2,sum)
PAIu2<-vegp$PAI[sel2,]
if (length(sel2) > 1) PAIu2<-apply(PAIu2,2,sum)
if (length(wgt2) > 1) {
PAIu<-PAIu1+(wgt1/(wgt1+wgt2))*(PAIu2-PAIu1)
} else PAIu<-PAIu1
} else {
zu<-z[length(z)]
wgt<-(hgt-reqhgt)/(hgt-zu)
PAIu<-wgt*vegp$PAI[length(z),]
}
dif<-abs(z-reqhgt)
sel<-which(dif==min(dif))[1]
leafdens<-vegp$PAI[sel,]/(z[2]-z[1])
} else PAIu<-rep(0,length(PAIt))
# Calculate wind speed 2 m above canopy
if (metopen) {
if (windhgt != 2) u <- u*4.87/log(67.8*windhgt-5.42)
u2<-u*log(67.8*(hgt+2)-5.42)/log(67.8*2-5.42)
} else {
u2 <- u
if (windhgt != (2+hgt)) u2 <- windprofile(u, windhgt, hgt+2, a = 2, PAIt, hgt)
}
u2[u2<0.5]<-0.5
cp<-cpair(tair)
ph<-phair(tair,pk)
zm<-ifelse(hgt>snowdep,roughlength(hgt, PAIt),0.003)
d<-ifelse(hgt>snowdep,zeroplanedis(hgt, PAIt),snowdep-0.003)
zh<-0.2*zm
hgt2<-ifelse(hgt>snowdep,hgt,snowdep)
a1<-attencoef(hgt,PAIt,vegp$cd,mean(vegp$iw))
a2<-attencoef(hgt,0,vegp$cd,mean(vegp$iw))
uz1<-.windprofile(u2,hgt+2,reqhgt,a1,hgt,PAIt)
uz2<-.windprofile(u2,snowdep+2,reqhgt,a2,0,0)
uz<-ifelse(hgt>snowdep,uz1,uz2)
uf<-(0.4*u2)/log((hgt2+2-d)/zm)
uf[uf<0.065]<-0.065
hgt2<-ifelse(hgt<snowdep,snowdep,hgt)
uh<-(uf/0.4)*log((hgt2-d)/zm)
uh[uh<uf]<-uf[uh<uf]
sas <- which(reqhgt>=snowdep)
sbs <- which(reqhgt<snowdep)
# Below snow
if (length(sbs)>0) {
tz2<-.snowtempf(snowdep,snow,reqhgt)[sbs]
rh2<-rep(100,length(tz2))
Rsw2<-rep(0,length(tz2))
Rlw2<-5.67*10^-8*0.85*(tz2+273.15)^4
tleaf2<-tz2[sbs]
}
# Calculate things that are common to below and above
leafdens<-(PAIt/hgt)*1.2
dp<-climdata$difrad/climdata$swrad
dp[is.na(dp)]<-0.5
dp[is.infinite(dp)]<-0.5
dp[dp>1]<-1
gtt<-gturb(u2,hgt+2,hgt+2,hgt,hgt,PAIt,tair,pk=pk)
# Above canopy
if (reqhgt >= hgt) {
# Calculate conductivities
gt0<-gcanopy(uh,hgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*2*0.71,uh,tair,5,pk,5)
gC<-layercond(climdata$swrad,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uh
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,snowtemp,lat,long,PAIt,0,pLAI,vegp$x,0.95,0.95,
0.95,0.85,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,snowtemp,relhum,pk,0.5,gtt,gt0,gha,gv,gL,Rabs,0.85,soilp$b,
soilp$psi_e,soilp$Smax,1,leafdens)
tn<-Th$tn
tleaf<-Th$tleaf
if (reqhgt == hgt) {
tz<-tn
rh<-Th$rh
} else {
# Calculate temperature and relative humidity at height z
gtc<-gturb(u2,hgt+2,reqhgt,hgt,hgt,PAIt,tair,pk=pk)
gt2<-1/(1/gtt-1/gtc)
eref<-(relhum/100)*satvap(tair)
ecan<-(Th$rh/100)*satvap(tn)
ea<-(gt2*eref+gtc*ecan)/(gt2+gtc)
tz<-(gt2*tair+gtc*tn)/(gt2+gtc)
rh<-(ea/satvap(tz))*100
}
rh[rh>100]<-100
Rsw<-climdata$swrad
sb<-5.67*10^-8
Rlw<-(1-climdata$skyem)*0.85*sb*(tz+273.15)^4
} else {
# Calculate conductivities
gtc<-gcanopy(uh,hgt,reqhgt,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gtt<-1/(1/gtt+1/gtc)
gt0<-gcanopy(uh,reqhgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*0.71*2,uz,tair,5,pk,5)
PAR<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=0.6)
gC<-layercond(PAR,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uz
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,snowtemp,lat,long,PAIt,PAIu,pLAI,vegp$x,0.95,0.95,
0.95,0.85,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,snowtemp,relhum,pk,0.5,gtt,gt0,gha,gv,gL,Rabs,vegp$vegem,soilp$b,
soilp$psi_e,soilp$Smax,1,leafdens)
tz<-Th$tn
tleaf<-Th$tleaf
rh<-Th$rh
# Radiation
Rsw<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=0.95)
# Calculate diffuse proportion
Rdif<-cantransdif(l=PAIu,ref=0.95)*climdata$swrad*dp
dp<-Rdif/Rsw
dp[is.na(dp)]<-1
dp[is.infinite(dp)]<-1
dp[dp<0]<-0
dp[dp>1]<-1
Rlw<-canlw(tair,PAIu,0.85,climdata$skyem,vegp$clump)$lwin
}
if (length(sbs)>0) {
tz[sbs]<-tz2
tleaf[sbs]<-tleaf2
rh[sbs]<-rh2
uz[sbs]<-0
Rsw[sbs]<-Rsw2
Rlw[sbs]<-Rlw2
}
metout<-data.frame(obs_time=climdata$obs_time,Tref=climdata$temp,Tloc=tz,tleaf=tleaf,
RHref=relhum,RHloc=rh,RSWloc=Rsw,RLWloc=Rlw,windspeed=uz,dp=dp)
return(metout)
}
#' Run model under steady state conditions
#' @description Rapid method for calculating below or above canopy or below ground microclimate
#' under steady-state at one user specified height.
#' @param climdata data.frame of climate variables needed to run the run the model (dataset should follow format of [weather()])
#' @param vegp a list of vegetation parameters as returned by [microctools::habitatvars()].
#' @param soilp a list of soil parameters as returned by [soilinit()]
#' @param nmrout a list of putputs from NicheMapR as returned by [runNMR()].
#' @param reqhgt height (m) for which microclimate is needed.
#' @param lat latitude of location (decimal degrees).
#' @param long longitude of location (decimal degrees).
#' @param metopen optional logical indicating whether the wind measurement used as an input to
#' the model is from a nearby weather station located in open ground (TRUE) or above the canopy
#' for which temperatures are modelled (FALSE - see details)
#' @param windhgt height above ground of wind measurement. If `metopen` is FALSE, must be above
#' canopy.
#' @param surfwet proportion of leaf surface acting as free water surface
#' @param groundem thermal emissivity of ground layer
#' @return a data.frame with the following columns:
#' @return `obs_time` time of observation. Same as in climdata.
#' @return `Tref` Temperature (deg C) at reference height as in climdata
#' @return `Tloc` Temperature (deg C) at height `reqhgt`
#' @return `tleaf` Leaf temperature (deg C) at height `reqhgt`. -999 if `reqhgt` above canopy
#' or below ground
#' @return `RHref` Relative humidity (percentage) at reference height as in climdata.
#' @return `RHloc` Relative humidity (percentage) at height `reqhgt`
#' @return `RSWloc` Total incoming shortwave radiation at height `reghgt` (W/m^2)
#' @return `RLWloc` Total downward longwave radiation at height `reqhgt` (W/m^2)
#' @return `windspeed` wind speed at height `reqhgt` (m/s)
#' @return `dp` fraction of Rswloc that is diffuse radiation
#'
#' @import microctools
#' @export
#'
#' @details This is a rapid implementation of model when time increments are hourly such that
#' transient heat fluxes and heat storage in canopy can be ignored. Computations are performed
#' simultaniously on all data, negating need to run in timesteps. Also includes implementation of
#' model with snow present. Should generally be run using wrapper function [runwithNMR()] but
#' provided as a standalone function in case data for multiple heights are needed, in whihc case
#' it can be run multiple times without also running NicheMapR.
runmodelS <- function(climdata, vegp, soilp, nmrout, reqhgt, lat, long, metopen = TRUE, windhgt = 2,
surfwet = 1, groundem = 0.95) {
if (!metopen & windhgt < vegp$hgt) {
stop("When metopen is FALSE, windhgt must be above canopy\n")
}
# (1) Unpack variables
tme<-as.POSIXlt(climdata$obs_time)
tair<-climdata$temp
relhum<-climdata$relhum
pk<-climdata$pres
u<-climdata$windspeed
hgt<-vegp$hgt
# Ensure at least some PAI
vegp$PAI[vegp$PAI<0.0001]<-0.0001
# Esimate total PAI and proportion LAI
PAIt<-apply(vegp$PAI,2,sum)
LAI<-(vegp$pLAI*vegp$PAI)
LAI<-apply(LAI,2,mean)
pLAI<-LAI/PAIt
m<-length(vegp$iw)
z<-c((1:m)-0.5)/m*hgt
# Esimate PAI above point and PAI of layer
if (reqhgt < hgt) {
sel<-which(z>reqhgt)
if (length(sel) > 1) {
wgt1<-abs(z[sel[1]]-reqhgt)
wgt2<-abs(z[sel[1]-1]-reqhgt)
sel2<-c(sel[1]-1,sel)
PAIu1<-vegp$PAI[sel,]
if (length(sel) > 1) PAIu1<-apply(PAIu1,2,sum)
PAIu2<-vegp$PAI[sel2,]
if (length(sel2) > 1) PAIu2<-apply(PAIu2,2,sum)
if (length(wgt2) > 1) {
PAIu<-PAIu1+(wgt1/(wgt1+wgt2))*(PAIu2-PAIu1)
} else PAIu<-PAIu1
} else {
zu<-z[length(z)]
wgt<-(hgt-reqhgt)/(hgt-zu)
PAIu<-wgt*vegp$PAI[length(z),]
}
dif<-abs(z-reqhgt)
sel<-which(dif==min(dif))[1]
leafdens<-vegp$PAI[sel,]/(z[2]-z[1])
} else PAIu<-rep(0,length(PAIt))
# Calculate wind speed 2 m above canopy
if (metopen) {
if (windhgt != 2) u <- u*4.87/log(67.8*windhgt-5.42)
u2<-u*log(67.8*(hgt+2)-5.42)/log(67.8*2-5.42)
} else {
u2 <- u
if (windhgt != (2+hgt)) u2 <- windprofile(u, windhgt, hgt+2, a = 2, PAIt, hgt)
}
# Get an approximate value for H for diabatic terms
cp<-cpair(tair)
ph<-phair(tair,pk)
u2[u2<0.5]<-0.5
zm<-roughlength(hgt, PAIt)
zh<-0.2*zm
d<-zeroplanedis(hgt, PAIt)
uf <- (0.4*u2)/log((2+hgt-d)/zm)
sb<-5.67*10^-8
Rlw<-(1-climdata$skyem)*sb*vegp$vegem*(tair+273.15)^4
Rnet<-(1-0.23)*climdata$swrad-Rlw
H<-0.2*Rnet
dba <- diabatic_cor(tair,pk,H,uf,hgt+2,d)
dba$psi_m<-ifelse(dba$psi_m>2.5,2.5,dba$psi_m)
dba$psi_h<-ifelse(dba$psi_h>2.5,2.5,dba$psi_h)
dba$psi_m<-ifelse(dba$psi_m< -2.5,-2.5,dba$psi_m)
dba$psi_h<-ifelse(dba$psi_h< -2.5,-2.5,dba$psi_h)
# Wind speed at user height
a<-attencoef(hgt,PAIt,vegp$cd,mean(vegp$iw))
uz<-.windprofile(u2,hgt+2,reqhgt,a,hgt,PAIt)
ln2 <- suppressWarnings(log((hgt - d) / zm))
ln2[ln2<0.55]<-0.55
uh <- (uf/0.4)*ln2
# Calculate things that are common to below and above
soilt<-nmrout$soiltemps
soilm<-nmrout$soilmoist
tground<- soilt$D0cm
theta<-soilm$WC0cm
leafdens<-PAIt/hgt
dp<-climdata$difrad/climdata$swrad
dp[is.na(dp)]<-0.5
dp[is.infinite(dp)]<-0.5
dp[dp>1]<-1
gtt<-gturb(u2,hgt+2,hgt+2,hgt,hgt,PAIt,tair,dba$psi_m,dba$psi_h,0.004,pk)
# Above canopy
if (reqhgt >= hgt) {
# Calculate conductivities
gt0<-gcanopy(uh,hgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*0.71,uh,tair,5,pk,5)
gC<-layercond(climdata$swrad,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uh
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,tground,lat,long,PAIt,0,pLAI,vegp$x,vegp$refls,vegp$refw,
vegp$refg,vegp$vegem,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,tground,relhum,pk,theta,gtt,gt0,gha,gv,gL,Rabs,vegp$vegem,soilp$b,
soilp$psi_e,soilp$Smax,surfwet,leafdens)
tn<-Th$tn
tleaf<-Th$tleaf
if (reqhgt == hgt) {
tz<-tn
rh<-Th$rh
} else {
# Calculate temperature and relative humidity at height z
gtc<-gturb(u2,hgt+2,reqhgt,hgt,hgt,PAIt,tair,dba$psi_m,dba$psi_h,0.004,pk)
gt2<-1/(1/gtt-1/gtc)
eref<-(relhum/100)*satvap(tair)
ecan<-(Th$rh/100)*satvap(tn)
ea<-(gt2*eref+gtc*ecan)/(gt2+gtc)
tz<-(gt2*tair+gtc*tn)/(gt2+gtc)
rh<-(ea/satvap(tz))*100
}
rh[rh>100]<-100
Rsw<-climdata$swrad
} else {
# Calculate conductivities
gtc<-gcanopy(uh,hgt,reqhgt,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gtt<-1/(1/gtt+1/gtc)
gt0<-gcanopy(uh,reqhgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*0.71,uz,tair,5,pk,5)
PAR<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=0.2)
gC<-layercond(PAR,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uz
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,tground,lat,long,PAIt,PAIu,pLAI,vegp$x,vegp$refls,vegp$refw,
vegp$refg,vegp$vegem,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,tground,relhum,pk,theta,gtt,gt0,gha,gv,gL,Rabs,vegp$vegem,soilp$b,
soilp$psi_e,soilp$Smax,surfwet,leafdens)
tz<-Th$tn
tleaf<-Th$tleaf
rh<-Th$rh
# Radiation
Rsw<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=vegp$refls)
# Calculate diffuse proportion
Rdif<-cantransdif(l=PAIu,ref=vegp$refls)*climdata$swrad*dp
dp<-Rdif/Rsw
dp[is.na(dp)]<-1
dp[is.infinite(dp)]<-1
dp[dp<0]<-0
dp[dp>1]<-1
Rlw<-canlw(tair,PAIu,1-vegp$vegem,climdata$skyem,vegp$clump)$lwin
}
metout<-data.frame(obs_time=climdata$obs_time,Tref=climdata$temp,Tloc=tz,tleaf=tleaf,
RHref=relhum,RHloc=rh,RSWloc=Rsw,RLWloc=Rlw,windspeed=uz,dp=dp)
# Consider snow
snow<-nmrout$snow
if (class(snow)[1] == "data.frame") {
snowtemp<-nmrout$SN1
} else snowtemp<-rep(0,length(tz))
mo<-nmrout$metout
snowdep<-mo$SNOWDEP
if (max(snowdep) > 0) {
mos <- .runmodelsnow(climdata,vegp,soilp,nmrout,reqhgt,lat,long,metopen,windhgt)
sel<-which(snowdep>0)
metout$Tloc[sel]<-mos$Tloc[sel]
metout$tleaf[sel]<-mos$tleaf[sel]
metout$RHloc[sel]<-mos$RHloc[sel]
metout$RSWloc[sel]<-mos$RSWloc[sel]
metout$RLWloc[sel]<-mos$RLWloc[sel]
metout$windspeed[sel]<-mos$windspeed[sel]
metout$dp[sel]<-mos$dp[sel]
}
return(metout)
}
#' Runs microclimate model in hourly timesteps with NicheMapR
#'
#' @description This function is a wrapper function for [runNMR()] and [runmodelS()] enabling
#' full microclimate model at any height above or below ground rapidly. NicheMapR is to
#' compute snowdepths, soil moistures and below-ground soil temperatures. The function
#' runmodelS is used to compute below or above canopy air temperatures.
#'
#' @param climdata data.frame of climate variables needed to run the run the model. The dataset should follow format
#' of [weather()]). Times must be in UTC.
#' @param prec vector of hourly or daily precipitation (mm) (see details).
#' @param vegp a list of vegetation parameters as returned by [microctools::habitatvars()].
#' @param soilp a list of soil parameters as returned by [soilinit()].
#' @param reqhgt height (m) for which microclimate is needed.
#' @param lat latitude of location (decimal degrees).
#' @param long longitude of location (decimal degrees).
#' @param altt elevation of location (m).
#' @param slope slope of ground surface (decimal degrees).
#' @param aspect aspect of ground surface (decimal degrees, 0 = north).
#' @param metopen optional logical indicating whether the wind measurement used as an input to
#' the model is from a nearby weather station located in open ground (TRUE) or above the canopy
#' for which temperatures are modelled (FALSE - see details)
#' @param windhgt height above ground of wind measurement. If `metopen` is FALSE, must be above
#' canopy.
#' @param surfwet proportion of leaf surface acting as free water surface.
#' @param groundem thermal emissivity of ground layer.
#' @param ERR Integrator error tolerance for soil temperature calculations.
#' @param cap Is organic cap present on soil surface? (1 = Yes, 0 = No).
#' @param hori Horizon angles (degrees), from 0 degrees azimuth (north) clockwise in 10 degree intervals.
#' @param maxpool Max depth for water pooling on the surface (mm), to account for runoff.
#' @param rainmult Rain multiplier for surface soil moisture (used to induce runon).
#' @param SoilMoist_Init Initial volumetric soil water content at each soil node (m^3/m^3)
#' @return a list of two objects:
#' @return a data.frame of microclimate variables for height `reqhgt` as returned by [runmodelS()]
#' @return a list of NicheMapR outputs as returned by [runNMR()]
#'
#' @import microctools
#' @import NicheMapR
#' @export
#' @seealso [runmodelS()], [runNMR()]
#'
#' @details This is a wrapper function for [runmodelS()] and [runNMR()] that allows fairly rapid
#' calaulation of microclimatic conditions below ground and/or above canopy. It uses the NicheMapR package
#' to calculate soil moisture-dependent soil temperatures, treating the vegetation canopy as a
#' single layer and then calculates air and leaf temperatures,air humidity and radiation using
#' runmodelS. Precipitation values are needed for computation of soil moisture and can either be
#' supplied as a vector of hourly values (must have identical number of values as in `climdata`) or
#' daily values (`climdata` must contain climate for entire days). It compares the length of the vector
#' to `climdata` to determine whether values are daily or hourly. In this function constant soil
#' poperties are are assumed throughout the soil profile. For more flexible implementation, [runNMR()]
#' and [runmodelS()] can be run seperately.
#'
#' @examples
#' require(NicheMapR)
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Run model using inbuilt weather datasets ~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' tme<-as.POSIXlt(weather$obs_time,tz="UTC")
#' vegp <- habitatvars(4, 50.2178, -5.32656, tme, m = 10) # Decidious broadleaf forest
#' soilp<- soilinit("Loam")
#' climrun <- runwithNMR(weather,dailyprecip,vegp,soilp,1,50.2178,-5.32656)
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Extract and view data from outputs ~~~~~~~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' metout<-climrun$metout
#' nmrout<-climrun$nmrout
#' soiltemps<-as.data.frame(nmrout$soil)
#' head(metout)
#' head(soiltemps)
#' head(nmrout$soilmoist)
#' head(nmrout$sunsnow)
#' head(nmrout$plant)
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Plot soil and air temperatures ~~~~~~~~~~~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' tsoil1<-soiltemps$D2.5cm
#' tsoil2<-soiltemps$D50cm
#' tair<-metout$Tloc
#' tref<-metout$Tref
#' tleaf<-metout$tleaf
#' Month<-as.POSIXct(metout$obs_time,tz="UTC")
#' tmn<-min(tsoil1,tsoil2,tair,tref,tleaf)
#' tmx<-max(tsoil1,tsoil2,tair,tref,tleaf)
#' par(mfrow=c(2,1))
#' # Air
#' plot(tref~Month,type="l",col=rgb(0,0,0,0.3),ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tleaf~Month,type="l",col=rgb(0,1,0,0.3),xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tair~Month,type="l",col=rgb(1,0,0,0.3),xlab="",ylab="",ylim=c(tmn,tmx))
#' # Soil
#' plot(tref~Month,type="l",col=rgb(0,0,0,0.5),ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil1~Month,type="l",col="red",xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil2~Month,type="l",col="blue",lwd=2,xlab="",ylab="",ylim=c(tmn,tmx))
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Plot data by monthly mean ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' tme<-as.POSIXlt(metout$obs_time,tz="UTC")
#' tsoil1m<-aggregate(tsoil1,by=list(tme$mon),mean)$x
#' tsoil2m<-aggregate(tsoil2,by=list(tme$mon),mean)$x
#' tairm<-aggregate(tair,by=list(tme$mon),mean)$x
#' trefm<-aggregate(tref,by=list(tme$mon),mean)$x
#' tleafm<-aggregate(tleaf,by=list(tme$mon),mean)$x
#' tmn<-min(tsoil1m,tsoil2m,tairm,trefm,tleafm)
#' tmx<-max(tsoil1m,tsoil2m,tairm,trefm,tleafm)
#' Month<-c(1:12)
#' # Air
#' par(mfrow=c(2,1))
#' plot(trefm~Month,type="l",lwd=2,col="darkgray",ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tleafm~Month,type="l",lwd=2,col="darkgreen",xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tairm~Month,type="l",lwd=2,col="red",xlab="",ylab="",ylim=c(tmn,tmx))
# Soil
#' plot(trefm~Month,type="l",lwd=2,col="darkgray",ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil1m~Month,type="l",lwd=2,col="red",xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil2m~Month,type="l",lwd=2,col="blue",xlab="",ylab="",ylim=c(tmn,tmx))
#' # ====================================================================== #
runwithNMR <- function(climdata, prec, vegp, soilp, reqhgt, lat, long, altt = 0, slope = 0,
aspect = 0, metopen = TRUE, windhgt = 2, surfwet = 1, groundem = 0.95,
ERR = 1.5, cap = 1, hori = rep(0,36), maxpool =1000, rainmult = 1,
SoilMoist_Init = c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3),
animal = FALSE) {
if (!metopen & windhgt < vegp$hgt) {
stop("When metopen is FALSE, windhgt must be above canopy\n")
}
# (1) Unpack variables
tair<-climdata$temp
relhum<-climdata$relhum
pk<-climdata$pres
u<-climdata$windspeed
hgt<-vegp$hgt
# (1) Run NicheMapR
PAIt<-apply(vegp$PAI,2,sum)
LAI<-(vegp$pLAI*vegp$PAI)
LAI<-apply(LAI,2,mean)
pLAI<-LAI/PAIt
LREFL<-mean(pLAI*vegp$refls+(1-pLAI)*vegp$reflp)
soiltype<-as.character(soilp$Soil.type)
DEP<-c(0,2.5,5,10,15,20,30,50,100,200)
PE = rep(1.1, 19)
KS = rep(0.0037, 19)
BB = rep(4.5, 19)
BD = rep(1.3, 19)
DD = rep(2.65, 19)
nmrout<-runNMR(climdata,prec,lat,long,1.95,hgt,2,PAIt,vegp$x,pLAI,
vegp$clump,vegp$refg,LREFL,groundem,DEP,altt,slope,aspect,
ERR,soiltype,PE,KS,BB,BD,DD,cap,hori,maxpool,rainmult,SoilMoist_Init,
animal=animal)
if (reqhgt < 0) {
dep<-c(0,2.5,5,10,15,20,30,50,100,200)
soilz<- -reqhgt*100
dif<-abs(soilz-dep)
o<-order(dif)
dif1<-dif[o[1]]
dif2<-dif[o[2]]
tz1<-soilt[,o[1]+1]
tz2<-soilt[,o[2]+1]
tz<-tz1*(dif2/(dif1+dif2))+tz2*(dif1/(dif1+dif2))
metout <- data.frame(obs_time=climdata$obs_time,Tref=tair,Tloc=tz,
tleaf=-999,RHref=relhum,RHloc=-999)
} else if (reqhgt < (hgt+2)) {
metout <- runmodelS(climdata,vegp,soilp,nmrout,reqhgt,lat,long,metopen,windhgt,surfwet,groundem)
} else {
metout <- data.frame(obs_time=climdata$obs_time,Tref=tair,Tloc=tair,
tleaf=-999,RHref=relhum,RHloc=relhum)
warning("Height out of range. Output climate identical to input")
}
climdata$temp <- metout$Tloc
climdata$relhum <- metout$RHloc
climdata$swrad <- metout$RSWloc
climdata$windspeed <- metout$windspeed
nmrout2<-runNMR(climdata=climdata,prec=prec,lat=lat,long=long,
Usrhyt=1.95,Veghyt=hgt,Refhyt=reqhgt,PAIt,vegp$x,pLAI,
vegp$clump,vegp$refg,LREFL,groundem,DEP,altt,slope,aspect,
ERR,soiltype,PE,KS,BB,BD,DD,cap,hori,maxpool,rainmult,SoilMoist_Init,
animal=animal)
nmrouts<-nmrout2$nmrout
nmrouts$soil<-nmrout$nmrout$soil
nmrouts$soilmoist<-nmrout$nmrout$moist
return(list(metout = metout, nmrout = nmrouts))
}
ilyamaclean/microclimc documentation built on July 28, 2023, 1:40 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions runwithNMR runmodelS .runmodelsnow .snowtempf tleafS .windprofile .leafabs2 runNMR
Documented in .leafabs2 runmodelS .runmodelsnow runNMR runwithNMR .snowtempf tleafS .windprofile
#' Runs NicheMapR model
#'
#' @description Wrapper function for running NicheMapR to derive soil moistures
#' @param climdata a data.frame of hourly weather variables in same format as [microclimc::weather()].
#' Column obs_time must give times in GMT/UTC.
#' @param prec a vector of daily or hourly precipitation values (mm).
#' @param lat latitude (decimal degrees, positive in northern hemisphere).
#' @param long longitude (decimal degrees, negative west of Greenwich meridian).
#' @param Usrhyt Local height (m) at which air temperature, wind speed and humidity are to be computed (see details).
#' @param Veghyt At of vegetation canopy (see details).
#' @param Refhyt Reference height (m) at which air input climate variables are measured.
#' @param PAI single value or vector of hourly or daily values of total plant area per unit ground area of canopy.
#' @param LOR Campbell leaf angle distrubution coefficient
#' @param pLAI fraction of PAI that is green vegetation.
#' @param clump clumpiness factor for canopy (0 to 1, 0 = even)
#' @param REFL A single numeric value of ground reflectivity of shortwave radiation.
#' @param LREFL A single numeric value of average canopy reflectivity of shortwave radiation.
#' @param SLE Thermal emissivity of ground
#' @param DEP Soil depths at which calculations are to be made (cm), must be 10 values
#' starting from 0, and more closely spaced near the surface.
#' @param ALTT elevation (m) of location.
#' @param SLOPE slope of location (decimal degrees).
#' @param ASPECT aspect of location (decimal degrees, 0 = north).
#' @param ERR Integrator error tolerance for soil temperature calculations.
#' @param soiltype soil type as for [microclimc::soilparams()]. Set to NA to use soil parameters
#' specified below.
#' @param PE Air entry potentials (J/kg) (19 values descending through soil for specified soil
#' nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param KS Saturated conductivities (kg s/m^3) (19 values descending through soil for
#' specified soil nodes in parameter DEP and points half way between). Ignored if soil type
#' is given.
#' @param BB Campbell's soil 'b' parameter (19 values descending through soil for specified
#' soil nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param BD Soil bulk density (Mg/m^3) (19 values descending through soil for specified
#' soil nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param DD Soil density (Mg/m^3) (19 values descending through soil for specified
#' soil nodes in parameter DEP and points half way between). Ignored if soil type is given.
#' @param cap Is organic cap present on soil surface? (1 = Yes, 0 = No).
#' @param hori Horizon angles (degrees), from 0 degrees azimuth (north) clockwise in 10 degree intervals.
#' @param maxpool Max depth for water pooling on the surface (mm), to account for runoff.
#' @param rainmult Rain multiplier for surface soil moisture (used to induce runon).
#' @param SoilMoist_Init Initial volumetric soil water content at each soil node (m^3/m^3)
#' @return a list with the following components:
#' @return `metout` a data.frame of microclimate variables for each hour at height `Usrhyt`.
#' @return `soiltemps` a data.frame of hourly soil temperatures for each depth node.
#' @return `soilmoist` a data.frame of hourly soil volumetric water content for each depth node.
#' @return `snowtemp` if snow present, a data,frame of snow temperature (°C), at each of the potential 8 snow layers (see details).
#' 0 if snow not present.
#' @return `plant` a data.frame of plant transpiration rates (g/m^2/hr), leaf water potentials
#' (J/kg) and root water potential (J/kg) at each of the 10 specified depths.
#' @return `nmrout` full NicheMapR output
#' @details Requires NicheMapR: devtools::install_github('mrke/NicheMapR'). NicheMapR is an integrated
#' soil moisture and temperature model that treats the vegetation as a single layer (https://mrke.github.io/).
#' This wrapper function, runs the NicheMapR::microclimate function with reduced parameter
#' inputs, by specifying sensible default values for other parameters where it is unlikely that these
#' would be known or explicitely measured at the site. The degree of canopy shading is worked out
#' explicitely from `PAI` at each hourly time interval by estimating canopy tranmission of direct
#' and diffuse radiation, and by adjusting the amount of radiation that would be absorbed
#' at the ground surface for a given slope and aspect. Roughness lengths and zero plane displacement
#' heights are limited to <2m in NicheMapR, so if `Veghyt` > 2, it is set to 2m.
#' This is adequate for computing soil moisture and total evapotranspiration, but will give false
#' air temperatures for heights below canopy (see [runwithNMR()]. If `soiltype` is given, the subsequent
#' soil parameters are computed from soil type. If `soiltype` is NA, soil paramaters must be specified.
#'
#' @author Ilya Maclean (i.m.d.maclean@exeter.ac.uk) and Urtzi Urzelai (urtzi.enriquez@gmail.com)
#' @import microctools
#' @import NicheMapR
#' @export
#' @examples
#' # Run NicheMapR with default parameters and inbuilt weather datasets
#' library(NicheMapR)
#' # Plot air temperatures
#' microout<-runNMR(weather,dailyprecip,50.2178,-5.32656,0.05,0.1,PAI=1)
#' metout <- microout$metout
#' tmn <- min(metout$TALOC,metout$TAREF)
#' tmx <- max(metout$TALOC,metout$TAREF)
#' dday <- metout$DOY+metout$TIME/1440 # Decimal day
#' plot(metout$TALOC~dday, type="l", col = "red", ylim=c(tmn,tmx), ylab = "Temperature")
#' par(new = T)
#' plot(metout$TAREF~dday, type="l", col = "blue", ylim=c(tmn,tmx), ylab = "", xlab = "")
#' # Plot soil temperatures
#' soiltemp<-microout$soiltemps
#' tmn <- min(soiltemp$D0cm,soiltemp$D30cm)
#' tmx <- max(soiltemp$D0cm,soiltemp$D30cm)
#' plot(soiltemp$D0cm~dday, type="l", col = "red", ylim=c(tmn,tmx), ylab = "Temperature")
#' par(new = T)
#' plot(soiltemp$D30cm~dday, type="l", col = "blue", ylim=c(tmn,tmx), ylab = "", xlab = "")
#' # Plot soil moistures
#' soilm <- microout$soilmoist
#' mmn <- min(soilm$WC2.5cm,soilm$WC30cm)
#' mmx <- max(soilm$WC2.5cm,soilm$WC30cm)
#' plot(soilm$WC2.5cm~dday, type="l", col = "red", ylim=c(mmn,mmx), ylab = "Soil moisture")
#' par(new = T)
#' plot(soilm$WC30cm~dday, type="l", col = "blue", ylim=c(mmn,mmx), ylab = "", xlab = "")
runNMR <- function(climdata, prec, lat, long, Usrhyt, Veghyt, Refhyt = 2, PAI = 3, LOR = 1,
pLAI = 0.8, clump = 0, REFL = 0.15, LREFL = 0.4, SLE = 0.95,
DEP = c(0,2.5,5,10,15,20,30,50,100,200), ALTT = 0, SLOPE = 0, ASPECT = 0,
ERR = 1.5, soiltype = "Loam", PE = NA, KS = NA, BB = NA, BD = NA, DD = NA,
cap = 1, hori = rep(0,36), maxpool = 1000, rainmult = 1,
SoilMoist_Init = c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3),
animal = FALSE) {
# Internal functions
.zeroplanedis<-function(h,pai) {
pai[pai<0.001]<-0.001
d<-(1-(1-exp(-sqrt(7.5*pai)))/sqrt(7.5*pai))*h
d
}
#' Calculate roughness length
.roughlength<-function(h,pai,d) {
Be<-sqrt(0.003+(0.2*pai)/2)
zm<-(h-d)*exp(-0.4/Be)
zm
}
# Hourly to daily
.dmaxmin <- function(x, fun) {
dx <- t(matrix(x, nrow = 24))
apply(dx, 1, fun)
}
# Get soil parameters
.getsoilparams<-function(soiltype) {
# Get soil parameters based on soil type
sel<-which(soilparams$Soil.type==soiltype)
BulkDensity<-soilparams$rho[sel]
SatWater <- soilparams$Smax[sel] # volumetric water content at saturation (0.1 bar matric potential) (m3/m3)
Clay <- CampNormTbl9_1$Clay[sel]*100 # clay con
KS<-rep(CampNormTbl9_1$Ks[sel],19)
BB<-rep(CampNormTbl9_1$b[sel],19)
PE<-rep(CampNormTbl9_1$airentry[sel],19)
return(list(BulkDensity=BulkDensity,SatWater=SatWater,Clay=Clay,KS=KS,BB=BB,PE=PE))
}
# Time and location parameters
ndays<-length(climdata$temp)/24
tme<-as.POSIXlt(climdata$obs_time,tz="UTC")
doy<-.dmaxmin(tme$yday+1,mean)
idayst <- 1 # start day (legacy parameter)
ida <- ndays # end day (legacy parameter)
HEMIS <- ifelse(lat < 0, 2, 1)
ALAT<-abs(trunc(lat))
AMINUT<-(abs(lat)-ALAT)*60
ALONG<- abs(trunc(long))
ALMINT<-(abs(long)-ALONG)*60
ALREF<-0
EC <- 0.0167238
# Air and wind vertical profile parameters
D0<-.zeroplanedis(Veghyt,mean(PAI))
RUF<-.roughlength(Veghyt,mean(PAI),D0)
D0<-0
ZH<-0.2*RUF
# Terrain and shading parameters
VIEWF<-1-sum(sin(hori*pi/180))/length(hori)
# Soil properties
Numtyps <- 2 # number of soil types
Nodes <- matrix(data = 0, nrow = 10, ncol = ndays) # array of all possible soil nodes
Nodes[1,1:ndays] <- 3 # deepest node for first substrate type
Nodes[2,1:ndays] <- 9 # deepest node for second substrate type
# Time varying environmental data
tannul <- mean(climdata$temp)
# Shade parameter
grasshade<-ifelse(Veghyt<0.5,1,0) # VT
if (length(PAI) == 1) {
PAI<-rep(PAI,ndays)
} else if (length(PAI) == ndays*24) {
PAI<-matrix(PAI,ncol=24,byrow=T)
PAI<-apply(PAI,1,mean)
}
if (length(PAI) != ndays) stop("PAI must be a single value or hourly/daily \n")
PAIc<-PAI^(1/(1-clump))
MINSHADES<-((1-exp(-PAIc))+clump)*100
MAXSHADES<-rep(100,ndays)
# Modal times of air temp, wind, humidity and cloud cover
TIMINS <- c(0, 0, 1, 1)
TIMAXS <- c(1, 1, 0, 0)
# Spinup
soilinit <- rep(tannul, 20)
spinup <- 1
# Decide whether to run snow model
if (length(prec) == ndays) {
RAINhr<-rep(prec/24,each=24)
} else RAINhr<-prec
snowmodel<-max(ifelse(climdata$temp<0,1,0)*ifelse(RAINhr>0,0,1))
densfun <- c(0.5979, 0.2178, 0.001, 0.0038)
intercept <- mean(MINSHADES) / 100 * 0.3
# Decide whether to run snowmodel
microinput<-c(ndays, RUF, ERR, Usrhyt, Refhyt, Numtyps, Z01 = 0, Z02 = 0, ZH1 = 0,
ZH2 = 0, idayst = 1, ida, HEMIS, ALAT, AMINUT, ALONG, ALMINT, ALREF,
SLOPE, ASPECT, ALTT, CMH2O = 1, microdaily = 1, tannul, EC, VIEWF,
snowtemp = 1.5, snowdens = 0.375, snowmelt = 0.9, undercatch = 1,
rainmult, runshade = 1, runmoist = 1, maxpool, evenrain = 1, snowmodel,
rainmelt = 0.0125, writecsv = 0, densfun, hourly = 1, rainhourly = 0,
lamb = 0, IUV = 0, RW = 2.5e+10, PC = -1500, RL = 2e+6, SP = 10, R1 = 0.001,
IM = 1e-06, MAXCOUNT = 500, IR = 0, message = 0, fail = 8760, snowcond = 0,
intercept, grasshade, solonly = 0, ZH, D0, TIMAXS, TIMINS, spinup,0,360)
tides <- matrix(data = 0, nrow = 24*ndays, ncol = 3)
SLES <- rep(SLE, ndays)
# Weather variables
TAIRhr<-climdata$temp
RHhr<-climdata$relhum
WNhr<-climdata$windspeed
PRESShr<-climdata$pres*1000
ea<-satvap(TAIRhr)*(RHhr/100)
eo<-1.24*(10*ea/(TAIRhr+273.15))^(1/7)
CLDhr<-((climdata$skyem-eo)/(1-eo))*100
CLDhr[CLDhr<0]<-0
CLDhr[CLDhr>100]<-100
TMAXX<-.dmaxmin(TAIRhr,max)
TMINN<-.dmaxmin(TAIRhr,min)
CCMAXX<-.dmaxmin(CLDhr,max)
CCMINN<-.dmaxmin(CLDhr,min)
RHMAXX<-.dmaxmin(RHhr,max)
RHMINN<-.dmaxmin(RHhr,min)
WNMAXX<-.dmaxmin(WNhr,max)
WNMINN<-.dmaxmin(WNhr,min)
PRESS<-.dmaxmin(PRESShr,min)
SOLRhr<-climdata$swrad*VIEWF
sb<-5.67*10^-8
IRDhr<-climdata$skyem*sb*(climdata$temp+273.15)^4
lt<-tme$hour+tme$min/60+tme$sec/3600
jd<-jday(tme=tme)
sa<-solalt(lt,lat,long,jd,0)
ZENhr<-90-sa
ZENhr[ZENhr>90]<-90
REFLS <- rep(REFL, ndays) # set up vector of soil reflectances for each day
# Soil surface wetness
RAIND<-.dmaxmin(RAINhr,max)
PCTWET<-ifelse(RAIND>0,90,0)
# Aerosol extinction coefficient profile
optdep.summer<-as.data.frame(rungads(lat,long,1,0))
optdep.winter<-as.data.frame(rungads(lat,long,1,0))
optdep<-cbind(optdep.winter[,1],rowMeans(cbind(optdep.summer[,2],optdep.winter[,2])))
optdep<-as.data.frame(optdep)
colnames(optdep)<-c("LAMBDA","OPTDEPTH")
a<-lm(OPTDEPTH~poly(LAMBDA,6,raw=TRUE),data=optdep)
LAMBDA<-c(290,295,300,305,310,315,320,330,340,350,360,370,380,390,400,420,440,460,480,
500,520,540,560,580,600,620,640,660,680,700,720,740,760,780,800,820,840,860,
880,900,920,940,960,980,1000,1020,1080,1100,1120,1140,1160,1180,1200,1220,
1240,1260,1280,1300,1320,1380,1400,1420,1440,1460,1480,1500,1540,1580,1600,
1620,1640,1660,1700,1720,1780,1800,1860,1900,1950,2000,2020,2050,2100,2120,
2150,2200,2260,2300,2320,2350,2380,2400,2420,2450,2490,2500,2600,2700,2800,
2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000)
TAI<-predict(a,data.frame(LAMBDA))
# Soil properties (mineral component only)
Thcond <- 2.5 # soil minerals thermal conductivity (W/mC)
SpecHeat <- 870 # soil minerals specific heat (J/kg-K)
if (class(DD) == "logical") {
Density<-2.560 # soil minerals density (Mg/m3)
} else Density<-mean(DD)
SP<-.getsoilparams(soiltype)
if (class(BD) == "logical") {
BulkDensity <- SP$BulkDensity # soil bulk density (kg/m3)
} else BulkDensity<-mean(BD)
SatWater<-SP$SatWater # volumetric water content at saturation (0.1 bar matric potential) (m3/m3)
Clay <- SP$Clay # clay content for matric potential calculations (%)
# Create matrix of soil properties
soilprops <- matrix(data = 0, nrow = 10, ncol = 6) # create an empty soil properties matrix
soilprops[1:2,1] <- BulkDensity # insert soil bulk density to profile 1
soilprops[1:2,2] <- SatWater # insert saturated water content to profile 1
soilprops[1:2,3] <- Thcond # insert thermal conductivity to profile 1
soilprops[1:2,4] <- SpecHeat # insert specific heat to profile 1
soilprops[1:2,5] <- Density # insert mineral density to profile 1
if(cap == 1) { # insert thermal conductivity to profile 1, and see if 'organic cap' added on top
soilprops[1,3]<-0.2 # mineral thermal conductivity
soilprops[1,4]<-1920
}
soilinit <- rep(tannul,20) # make inital soil temps equal to mean annual
if (class(PE) == "logical") PE<-SP$PE
if (class(KS) == "logical") KS<-SP$KS
if (class(BB) == "logical") BB<-SP$BB
if (class(BD) == "logical") BD<-rep(BulkDensity,19)
if (class(DD) == "logical") DD<-rep(Density,19)
# Initial soil moistures
moists <- matrix(nrow=10, ncol = ndays, data = 0) # set up an empty vector for soil moisture values through time
moists[1:10,] <- SoilMoist_Init # insert inital soil moisture
# Calculate monthly mean rainfall
RAINFALL<-.dmaxmin(RAINhr, sum)
tannulrun <- rep(tannul,ndays) # deep soil temperature (2m) (deg C)
# Root density at each node:
L<-c(0,0,8.2,8,7.8,7.4,7.1,6.4,5.8,4.8,4,1.8,0.9,0.6,0.8,0.4,0.4,0,0)*10000
LAI<-pLAI*PAI
micro<-list(microinput = microinput, tides=tides, doy = doy, SLES = SLES, DEP = DEP,
Nodes=Nodes, MAXSHADES = MAXSHADES, MINSHADES = MINSHADES, TMAXX = TMAXX,
TMINN = TMINN, RHMAXX = RHMAXX, RHMINN = RHMINN, CCMAXX = CCMAXX, CCMINN = CCMINN,
WNMAXX = WNMAXX, WNMINN = WNMINN, TAIRhr = TAIRhr, RHhr = RHhr, WNhr = WNhr,
CLDhr = CLDhr, SOLRhr = SOLRhr, RAINhr = RAINhr, ZENhr = ZENhr, IRDhr = IRDhr,
REFLS = REFLS, PCTWET = PCTWET, soilinit = soilinit, hori = hori,
TAI = TAI, soilprops = soilprops, moists = moists, RAINFALL = RAINFALL,
tannulrun = tannulrun, PE = PE, KS = KS, BB = BB, BD = BD, DD = DD,
L = L, LAI = LAI)
# Run model
microut<-microclimate(micro)
if(animal==FALSE){
metout<-as.data.frame(microut$metout)
soil<-as.data.frame(microut$soil)
soilmoist<-as.data.frame(microut$soilmoist)
plant<-as.data.frame(microut$plant)
if (snowmodel == 1) {
snow <- as.data.frame(microut$sunsnow)
} else snow <- 0
return(list(metout=metout,soiltemps=soil,soilmoist=soilmoist,snowtemp=snow,plant=plant,nmrout=microut))
}
if(animal == TRUE){
metout<-microut$metout
shadmet<-microut$shadmet
soil<-microut$soil
shadsoil<-microut$shadsoil
soilmoist<-microut$soilmoist
shadmoist <- microut$shadmoist
humid <- microut$humid
shadhumid <- microut$shadhumid
soilpot <- microut$soilpot
shadpot <- microut$shadpot
plant <- microut$plant
shadplant <- microut$shadplant
if (snowmodel == 1) {
snow <- microut$sunsnow
shdsnow <- microut$shdsnow
} else snow <- 0
drlam <- microut$drlam
drrlam <- microut$drrlam
srlam <- microut$srlam
timeinterval = 365
days <- rep(seq(1, timeinterval * nyears), 24)
days <- days[order(days)]
dates <- days + metout[, 2]/60/24 - 1
dates2 <- seq(1, timeinterval * nyears)
if (snowmodel == 1) {
nmrout_full <- list(soil = soil, shadsoil = shadsoil,
metout = metout, shadmet = shadmet, soilmoist = soilmoist,
shadmoist = shadmoist, humid = humid, shadhumid = shadhumid,
soilpot = soilpot, shadpot = shadpot, sunsnow = snow,
shdsnow = shdsnow, plant = plant, shadplant = shadplant,
RAINFALL = RAINFALL, ndays = ndays, elev = ALTT,
REFL = REFL[1], longlat = c(x[1], x[2]), nyears = nyears,
timeinterval = timeinterval, minshade = MINSHADES,
maxshade = MAXSHADES, DEP = DEP, drlam = drlam,
drrlam = drrlam, srlam = srlam, dates = dates,
dates2 = dates2)
}
else {
nmrout_full <- list(soil = soil, shadsoil = shadsoil,
metout = metout, shadmet = shadmet, soilmoist = soilmoist,
shadmoist = shadmoist, humid = humid, shadhumid = shadhumid,
soilpot = soilpot, shadpot = shadpot, plant = plant,
shadplant = shadplant, RAINFALL = RAINFALL,
ndays = ndays, elev = ALTT, REFL = REFL[1],
longlat = c(x[1], x[2]), nyears = nyears, timeinterval = timeinterval,
minshade = MINSHADES, maxshade = MAXSHADES,
DEP = DEP, drlam = drlam, drrlam = drrlam,
srlam = srlam, dates = dates, dates2 = dates2)
}
metout <- data.frame(metout)
soil <- data.frame(soil)
soilmoist <- data.frame(soilmoist)
snow <- data.frame(snow)
plant <- data.frame(plant)
return(list(metout=metout,soiltemps=soil,soilmoist=soilmoist,snowtemp=snow,plant=plant,nmrout=nmrout_full))
}
}
#' Internal function for calculating lead absorbed radiation on vector
.leafabs2 <-function(Rsw, tme, tair, tground, lat, long, PAIt, PAIu, pLAI, x, refls, refw, refg, vegem, skyem, dp,
merid = round(long/15, 0) * 15, dst = 0, clump = 0) {
jd<-jday(tme = tme)
lt<-tme$hour+tme$min/60+tme$sec/3600
sa<-solalt(lt,lat,long,jd,merid=merid,dst=dst)
sa2<-ifelse(sa<5,5,sa)
ref<-pLAI*refls+(1-pLAI)*refw
sunl<-psunlit(PAIu,x,sa,clump)
mul<-radmult(x,sa2)
aRsw <- (1-ref) * cansw(Rsw,dp,jd,lt,lat,long,PAIu,x,ref,merid=merid,dst=dst,clump=clump)
aGround <- refg * cansw(Rsw,dp,jd,lt,lat,long,PAIt,x,ref,merid=merid,dst=dst,clump=clump)
aGround <- (1-ref) * cansw(aGround,dp,jd,lt,lat,long,PAIt,x,ref,merid=merid,dst=dst,clump=clump)
aRsw <- sunl * mul * aRsw + (1-sunl) * aRsw + aGround
aRlw <- canlw(tair, PAIu, 1-vegem, skyem = skyem, clump = clump)$lwabs
return(aRsw+aRlw)
}
#' Internal function for calculating lead absorbed radiation on vector
.windprofile <- function(ui, zi, zo, a = 2, hgt, PAI, psi_m = 0, hgtg = 0.05 * hgt, zm0 = 0.004) {
d<-zeroplanedis(hgt,PAI)
zm<-roughlength(hgt,PAI,zm0)
ln2<-suppressWarnings(log((zi-d)/zm)+psi_m)
ln2[ln2<0.55]<-0.55
uf<-(0.4*ui)/ln2
if (zo >= hgt) {
ln2<-suppressWarnings(log((zo-d)/zm)+psi_m)
ln2[ln2<0.55]<-0.55
uo<-(uf/0.4)*ln2
} else {
ln2<-suppressWarnings(log((hgt-d)/zm)+psi_m)
ln2[ln2<0.55]<-0.55
uh<-(uf/0.4)*ln2
if (zi >= (0.1*hgt)) {
uo<-uh*exp(a*(zo/hgt)-1)
} else {
uo<-uh*exp(a*((0.1*hgt)/hgt)-1)
zmg<-0.1*hgtg
ln1<-log((0.1*hgt)/zmg)+psi_m
uf<-0.4*(uo/ln1)
ln2<-log(zo/zmg)+psi_m
uo<-(uf/0.4)*ln2
}
}
uo
}
#' Steady-state leaf and air temperature
#' @description Rapid method for calculating steady-state leaf and air temperature at one height
#' @param tair air temperature at reference height (deg C)
#' @param tground ground surface temperature (as returned by e.g. [runNMR()])
#' @param relhum relative humidity at reference height (percentage)
#' @param pk atmospheric pressure (kPa)
#' @param theta volumetric water content of upper most soil layer in current time step (m^3 / m^3)
#' as returned by e.g. [runNMR()]
#' @param gtt molar conductivity from leaf height to reference height as returned by [gturb()]
#' and [gcanopy()]
#' @param gt0 molar conductivity from leaf height to ground as returned by [gcanopy()]
#' @param gha boundary layer conductance of leaf as returned by [gforcedfree()]
#' @param gL combined boundary layer and leaf-air conductance given by 1/(1/gha+1/(uz*ph))
#' @param gv combined boundary layer and stomatal conductance of leaf
#' @param Rabs radiation absorbed by leaf as returned by e.g. [leafabs()]
#' @param soilb Shape parameter for Campbell soil model (dimensionless, > 1) as returned by
#' [soilinit()]
#' @param Psie Soil matric potential (J / m^3) as returned by [soilinit()]
#' @param Smax Volumetric water content at saturation (m3 / m3)
#' @param surfwet proportion of leaf surface acting as free water surface
#' @param leafdens Total one sided leaf area per m^3 at desired height
#' @return a list of the following:
#' @return `tleaf` leaf temperature
#' @return `tn` air temperature
#' @return `rh` relative humidity
#' @import microctools
#' @export
tleafS <- function(tair, tground, relhum, pk, theta, gtt, gt0, gha, gv, gL, Rabs, vegem, soilb,
Psie, Smax, surfwet, leafdens) {
###
cp<-cpair(tair)
# Air temperature expressed as leaf temperature
aL<-(gtt*(tair+273.15)+gt0*(tground+273.15))/(gtt+gt0)
bL<-(leafdens*gL)/(gtt+gt0)
# Vapour pressures
es<-satvap(tair)
eref <- (relhum/100)*es
rhsoil<-soilrh(theta,soilb,Psie,Smax,tground)
esoil<-rhsoil*satvap(tground)
# add a small correction to improve delta estimate
sb<-5.67*10^-8
Rnet<-Rabs-0.97*sb*(tair+273.15)^4
tle<-tair+(0.5*Rnet)/(cp*gha)
wgt1<-ifelse(leafdens<1,leafdens/2,1-0.5/leafdens)
wgt2<-1-wgt1
tapprox<-wgt1*tle+wgt2*tair
delta <- 4098*(0.6108*exp(17.27*tapprox/(tapprox+237.3)))/(tapprox+237.3)^2
ae<-(gtt*eref+gt0*esoil+gv*es)/(gtt+gt0+gv)
be<-(gv*delta)/(gtt+gt0+gv)
# Sensible heat
bH<-gha*cp
# Latent heat
lambda <- (-42.575*tair+44994)
aX<-((lambda*gv)/pk)*(surfwet*es-ae)
bX<-((lambda*gv)/pk)*(surfwet*delta-be)
aX[aX<0]<-0
bX[bX<0]<-0
# Emmited radiation
aR<-sb*vegem*aL^4
bR<-4*vegem*sb*(aL^3*bL+(tair+273.15)^3)
# Leaf temperature
dTL <- (Rabs-aR-aX)/(1+bR+bX+bH)
# tz pass 1
tn<-aL-273.15+bL*dTL
tleaf<-tn+dTL
# new vapour pressure
eanew<-ae+be*dTL
eanew[eanew<0.01]<-0.01
tmin<-dewpoint(eanew,tn,ice = TRUE)
esnew<-satvap(tn)
eanew<-ifelse(eanew>esnew,esnew,eanew)
rh<-(eanew/esnew)*100
# Set both tair and tleaf so as not to drop below dewpoint
tleaf<-ifelse(tleaf<tmin,tmin,tleaf)
tn<-ifelse(tn<tmin,tmin,tn)
tmax<-ifelse(tn+20<80,tn+20,80)
tleaf<-ifelse(tleaf>tmax,tmax,tleaf)
# cap upper limits of both tair and tleaf
tmx<-pmax(tground+5,tair+5)
tn<-ifelse(tn>tmx,tmx,tn)
tmx<-pmax(tground+20,tair+20)
tleaf<-ifelse(tleaf>tmx,tmx,tleaf)
# cap lower limits
tmn<-tair-7
tn<-ifelse(tn<tmn,tmn,tn)
tmn<-tair-20
tleaf<-ifelse(tleaf<tmn,tmn,tleaf)
return(list(tleaf=tleaf,tn=tn,rh=rh))
}
#' Internation function for computing snow temperature
.snowtempf <- function(snowdepth,snow,reqhgt) {
ndepths<-rev(c(0,2.5,5,10,20,50,100,200,300))/100
ha<-ndepths-reqhgt
selb<-which(ha<=0)[1]# layer below
sela<-selb-1 # layer below or equal
sm<-abs(reqhgt-ndepths[sela])+abs(reqhgt-ndepths[selb])
wgt1<-1-abs(reqhgt-ndepths[sela])/sm # weight above
wgt2<-1-abs(reqhgt-ndepths[selb])/sm # weight below
hs<-which(ndepths[sela]>(snowdepth/100))
wgt1<-rep(wgt1,length(snowdepth))
wgt2<-rep(wgt2,length(snowdepth))
wgt1[hs]<-0
wgt2[hs]<-1
stemp<-wgt1*snow[,sela+2]+wgt2*snow[,selb+2]
sel<-which(snowdepth<reqhgt)
stemp[sel]<-0
stemp
}
#' Internal function for running model with snow
.runmodelsnow <- function(climdata, vegp, soilp, nmrout, reqhgt, lat, long, metopen = TRUE, windhgt = 2) {
# Snow
snow<-nmrout$snow
if (class(snow)[1] == "data.frame") {
snowtemp<-snow$SN1
} else snowtemp<-rep(0,length(L))
metout<-nmrout$metout
snowdep<-metout$SNOWDEP/100
# (1) Unpack variables
tme<-as.POSIXlt(climdata$obs_time)
tair<-climdata$temp
relhum<-climdata$relhum
pk<-climdata$pres
u<-climdata$windspeed
hgt<-vegp$hgt
# Esimate PAIu
PAIt<-apply(vegp$PAI,2,sum)
PAIt[PAIt<0.001]<-0.001
LAI<-(vegp$pLAI*vegp$PAI)
LAI<-apply(LAI,2,mean)
pLAI<-LAI/PAIt
pLAI[is.na(pLAI)]<-0.8
m<-length(vegp$iw)
z<-c((1:m)-0.5)/m*hgt
# Esimate PAI above point and PAI of layer
if (reqhgt < hgt) {
sel<-which(z>reqhgt)
if (length(sel) > 1) {
wgt1<-abs(z[sel[1]]-reqhgt)
wgt2<-abs(z[sel[1]-1]-reqhgt)
sel2<-c(sel[1]-1,sel)
PAIu1<-vegp$PAI[sel,]
if (length(sel) > 1) PAIu1<-apply(PAIu1,2,sum)
PAIu2<-vegp$PAI[sel2,]
if (length(sel2) > 1) PAIu2<-apply(PAIu2,2,sum)
if (length(wgt2) > 1) {
PAIu<-PAIu1+(wgt1/(wgt1+wgt2))*(PAIu2-PAIu1)
} else PAIu<-PAIu1
} else {
zu<-z[length(z)]
wgt<-(hgt-reqhgt)/(hgt-zu)
PAIu<-wgt*vegp$PAI[length(z),]
}
dif<-abs(z-reqhgt)
sel<-which(dif==min(dif))[1]
leafdens<-vegp$PAI[sel,]/(z[2]-z[1])
} else PAIu<-rep(0,length(PAIt))
# Calculate wind speed 2 m above canopy
if (metopen) {
if (windhgt != 2) u <- u*4.87/log(67.8*windhgt-5.42)
u2<-u*log(67.8*(hgt+2)-5.42)/log(67.8*2-5.42)
} else {
u2 <- u
if (windhgt != (2+hgt)) u2 <- windprofile(u, windhgt, hgt+2, a = 2, PAIt, hgt)
}
u2[u2<0.5]<-0.5
cp<-cpair(tair)
ph<-phair(tair,pk)
zm<-ifelse(hgt>snowdep,roughlength(hgt, PAIt),0.003)
d<-ifelse(hgt>snowdep,zeroplanedis(hgt, PAIt),snowdep-0.003)
zh<-0.2*zm
hgt2<-ifelse(hgt>snowdep,hgt,snowdep)
a1<-attencoef(hgt,PAIt,vegp$cd,mean(vegp$iw))
a2<-attencoef(hgt,0,vegp$cd,mean(vegp$iw))
uz1<-.windprofile(u2,hgt+2,reqhgt,a1,hgt,PAIt)
uz2<-.windprofile(u2,snowdep+2,reqhgt,a2,0,0)
uz<-ifelse(hgt>snowdep,uz1,uz2)
uf<-(0.4*u2)/log((hgt2+2-d)/zm)
uf[uf<0.065]<-0.065
hgt2<-ifelse(hgt<snowdep,snowdep,hgt)
uh<-(uf/0.4)*log((hgt2-d)/zm)
uh[uh<uf]<-uf[uh<uf]
sas <- which(reqhgt>=snowdep)
sbs <- which(reqhgt<snowdep)
# Below snow
if (length(sbs)>0) {
tz2<-.snowtempf(snowdep,snow,reqhgt)[sbs]
rh2<-rep(100,length(tz2))
Rsw2<-rep(0,length(tz2))
Rlw2<-5.67*10^-8*0.85*(tz2+273.15)^4
tleaf2<-tz2[sbs]
}
# Calculate things that are common to below and above
leafdens<-(PAIt/hgt)*1.2
dp<-climdata$difrad/climdata$swrad
dp[is.na(dp)]<-0.5
dp[is.infinite(dp)]<-0.5
dp[dp>1]<-1
gtt<-gturb(u2,hgt+2,hgt+2,hgt,hgt,PAIt,tair,pk=pk)
# Above canopy
if (reqhgt >= hgt) {
# Calculate conductivities
gt0<-gcanopy(uh,hgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*2*0.71,uh,tair,5,pk,5)
gC<-layercond(climdata$swrad,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uh
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,snowtemp,lat,long,PAIt,0,pLAI,vegp$x,0.95,0.95,
0.95,0.85,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,snowtemp,relhum,pk,0.5,gtt,gt0,gha,gv,gL,Rabs,0.85,soilp$b,
soilp$psi_e,soilp$Smax,1,leafdens)
tn<-Th$tn
tleaf<-Th$tleaf
if (reqhgt == hgt) {
tz<-tn
rh<-Th$rh
} else {
# Calculate temperature and relative humidity at height z
gtc<-gturb(u2,hgt+2,reqhgt,hgt,hgt,PAIt,tair,pk=pk)
gt2<-1/(1/gtt-1/gtc)
eref<-(relhum/100)*satvap(tair)
ecan<-(Th$rh/100)*satvap(tn)
ea<-(gt2*eref+gtc*ecan)/(gt2+gtc)
tz<-(gt2*tair+gtc*tn)/(gt2+gtc)
rh<-(ea/satvap(tz))*100
}
rh[rh>100]<-100
Rsw<-climdata$swrad
sb<-5.67*10^-8
Rlw<-(1-climdata$skyem)*0.85*sb*(tz+273.15)^4
} else {
# Calculate conductivities
gtc<-gcanopy(uh,hgt,reqhgt,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gtt<-1/(1/gtt+1/gtc)
gt0<-gcanopy(uh,reqhgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*0.71*2,uz,tair,5,pk,5)
PAR<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=0.6)
gC<-layercond(PAR,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uz
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,snowtemp,lat,long,PAIt,PAIu,pLAI,vegp$x,0.95,0.95,
0.95,0.85,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,snowtemp,relhum,pk,0.5,gtt,gt0,gha,gv,gL,Rabs,vegp$vegem,soilp$b,
soilp$psi_e,soilp$Smax,1,leafdens)
tz<-Th$tn
tleaf<-Th$tleaf
rh<-Th$rh
# Radiation
Rsw<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=0.95)
# Calculate diffuse proportion
Rdif<-cantransdif(l=PAIu,ref=0.95)*climdata$swrad*dp
dp<-Rdif/Rsw
dp[is.na(dp)]<-1
dp[is.infinite(dp)]<-1
dp[dp<0]<-0
dp[dp>1]<-1
Rlw<-canlw(tair,PAIu,0.85,climdata$skyem,vegp$clump)$lwin
}
if (length(sbs)>0) {
tz[sbs]<-tz2
tleaf[sbs]<-tleaf2
rh[sbs]<-rh2
uz[sbs]<-0
Rsw[sbs]<-Rsw2
Rlw[sbs]<-Rlw2
}
metout<-data.frame(obs_time=climdata$obs_time,Tref=climdata$temp,Tloc=tz,tleaf=tleaf,
RHref=relhum,RHloc=rh,RSWloc=Rsw,RLWloc=Rlw,windspeed=uz,dp=dp)
return(metout)
}
#' Run model under steady state conditions
#' @description Rapid method for calculating below or above canopy or below ground microclimate
#' under steady-state at one user specified height.
#' @param climdata data.frame of climate variables needed to run the run the model (dataset should follow format of [weather()])
#' @param vegp a list of vegetation parameters as returned by [microctools::habitatvars()].
#' @param soilp a list of soil parameters as returned by [soilinit()]
#' @param nmrout a list of putputs from NicheMapR as returned by [runNMR()].
#' @param reqhgt height (m) for which microclimate is needed.
#' @param lat latitude of location (decimal degrees).
#' @param long longitude of location (decimal degrees).
#' @param metopen optional logical indicating whether the wind measurement used as an input to
#' the model is from a nearby weather station located in open ground (TRUE) or above the canopy
#' for which temperatures are modelled (FALSE - see details)
#' @param windhgt height above ground of wind measurement. If `metopen` is FALSE, must be above
#' canopy.
#' @param surfwet proportion of leaf surface acting as free water surface
#' @param groundem thermal emissivity of ground layer
#' @return a data.frame with the following columns:
#' @return `obs_time` time of observation. Same as in climdata.
#' @return `Tref` Temperature (deg C) at reference height as in climdata
#' @return `Tloc` Temperature (deg C) at height `reqhgt`
#' @return `tleaf` Leaf temperature (deg C) at height `reqhgt`. -999 if `reqhgt` above canopy
#' or below ground
#' @return `RHref` Relative humidity (percentage) at reference height as in climdata.
#' @return `RHloc` Relative humidity (percentage) at height `reqhgt`
#' @return `RSWloc` Total incoming shortwave radiation at height `reghgt` (W/m^2)
#' @return `RLWloc` Total downward longwave radiation at height `reqhgt` (W/m^2)
#' @return `windspeed` wind speed at height `reqhgt` (m/s)
#' @return `dp` fraction of Rswloc that is diffuse radiation
#'
#' @import microctools
#' @export
#'
#' @details This is a rapid implementation of model when time increments are hourly such that
#' transient heat fluxes and heat storage in canopy can be ignored. Computations are performed
#' simultaniously on all data, negating need to run in timesteps. Also includes implementation of
#' model with snow present. Should generally be run using wrapper function [runwithNMR()] but
#' provided as a standalone function in case data for multiple heights are needed, in whihc case
#' it can be run multiple times without also running NicheMapR.
runmodelS <- function(climdata, vegp, soilp, nmrout, reqhgt, lat, long, metopen = TRUE, windhgt = 2,
surfwet = 1, groundem = 0.95) {
if (!metopen & windhgt < vegp$hgt) {
stop("When metopen is FALSE, windhgt must be above canopy\n")
}
# (1) Unpack variables
tme<-as.POSIXlt(climdata$obs_time)
tair<-climdata$temp
relhum<-climdata$relhum
pk<-climdata$pres
u<-climdata$windspeed
hgt<-vegp$hgt
# Ensure at least some PAI
vegp$PAI[vegp$PAI<0.0001]<-0.0001
# Esimate total PAI and proportion LAI
PAIt<-apply(vegp$PAI,2,sum)
LAI<-(vegp$pLAI*vegp$PAI)
LAI<-apply(LAI,2,mean)
pLAI<-LAI/PAIt
m<-length(vegp$iw)
z<-c((1:m)-0.5)/m*hgt
# Esimate PAI above point and PAI of layer
if (reqhgt < hgt) {
sel<-which(z>reqhgt)
if (length(sel) > 1) {
wgt1<-abs(z[sel[1]]-reqhgt)
wgt2<-abs(z[sel[1]-1]-reqhgt)
sel2<-c(sel[1]-1,sel)
PAIu1<-vegp$PAI[sel,]
if (length(sel) > 1) PAIu1<-apply(PAIu1,2,sum)
PAIu2<-vegp$PAI[sel2,]
if (length(sel2) > 1) PAIu2<-apply(PAIu2,2,sum)
if (length(wgt2) > 1) {
PAIu<-PAIu1+(wgt1/(wgt1+wgt2))*(PAIu2-PAIu1)
} else PAIu<-PAIu1
} else {
zu<-z[length(z)]
wgt<-(hgt-reqhgt)/(hgt-zu)
PAIu<-wgt*vegp$PAI[length(z),]
}
dif<-abs(z-reqhgt)
sel<-which(dif==min(dif))[1]
leafdens<-vegp$PAI[sel,]/(z[2]-z[1])
} else PAIu<-rep(0,length(PAIt))
# Calculate wind speed 2 m above canopy
if (metopen) {
if (windhgt != 2) u <- u*4.87/log(67.8*windhgt-5.42)
u2<-u*log(67.8*(hgt+2)-5.42)/log(67.8*2-5.42)
} else {
u2 <- u
if (windhgt != (2+hgt)) u2 <- windprofile(u, windhgt, hgt+2, a = 2, PAIt, hgt)
}
# Get an approximate value for H for diabatic terms
cp<-cpair(tair)
ph<-phair(tair,pk)
u2[u2<0.5]<-0.5
zm<-roughlength(hgt, PAIt)
zh<-0.2*zm
d<-zeroplanedis(hgt, PAIt)
uf <- (0.4*u2)/log((2+hgt-d)/zm)
sb<-5.67*10^-8
Rlw<-(1-climdata$skyem)*sb*vegp$vegem*(tair+273.15)^4
Rnet<-(1-0.23)*climdata$swrad-Rlw
H<-0.2*Rnet
dba <- diabatic_cor(tair,pk,H,uf,hgt+2,d)
dba$psi_m<-ifelse(dba$psi_m>2.5,2.5,dba$psi_m)
dba$psi_h<-ifelse(dba$psi_h>2.5,2.5,dba$psi_h)
dba$psi_m<-ifelse(dba$psi_m< -2.5,-2.5,dba$psi_m)
dba$psi_h<-ifelse(dba$psi_h< -2.5,-2.5,dba$psi_h)
# Wind speed at user height
a<-attencoef(hgt,PAIt,vegp$cd,mean(vegp$iw))
uz<-.windprofile(u2,hgt+2,reqhgt,a,hgt,PAIt)
ln2 <- suppressWarnings(log((hgt - d) / zm))
ln2[ln2<0.55]<-0.55
uh <- (uf/0.4)*ln2
# Calculate things that are common to below and above
soilt<-nmrout$soiltemps
soilm<-nmrout$soilmoist
tground<- soilt$D0cm
theta<-soilm$WC0cm
leafdens<-PAIt/hgt
dp<-climdata$difrad/climdata$swrad
dp[is.na(dp)]<-0.5
dp[is.infinite(dp)]<-0.5
dp[dp>1]<-1
gtt<-gturb(u2,hgt+2,hgt+2,hgt,hgt,PAIt,tair,dba$psi_m,dba$psi_h,0.004,pk)
# Above canopy
if (reqhgt >= hgt) {
# Calculate conductivities
gt0<-gcanopy(uh,hgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*0.71,uh,tair,5,pk,5)
gC<-layercond(climdata$swrad,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uh
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,tground,lat,long,PAIt,0,pLAI,vegp$x,vegp$refls,vegp$refw,
vegp$refg,vegp$vegem,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,tground,relhum,pk,theta,gtt,gt0,gha,gv,gL,Rabs,vegp$vegem,soilp$b,
soilp$psi_e,soilp$Smax,surfwet,leafdens)
tn<-Th$tn
tleaf<-Th$tleaf
if (reqhgt == hgt) {
tz<-tn
rh<-Th$rh
} else {
# Calculate temperature and relative humidity at height z
gtc<-gturb(u2,hgt+2,reqhgt,hgt,hgt,PAIt,tair,dba$psi_m,dba$psi_h,0.004,pk)
gt2<-1/(1/gtt-1/gtc)
eref<-(relhum/100)*satvap(tair)
ecan<-(Th$rh/100)*satvap(tn)
ea<-(gt2*eref+gtc*ecan)/(gt2+gtc)
tz<-(gt2*tair+gtc*tn)/(gt2+gtc)
rh<-(ea/satvap(tz))*100
}
rh[rh>100]<-100
Rsw<-climdata$swrad
} else {
# Calculate conductivities
gtc<-gcanopy(uh,hgt,reqhgt,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gtt<-1/(1/gtt+1/gtc)
gt0<-gcanopy(uh,reqhgt,0,tair,tair,hgt,PAIt,vegp$cd,mean(vegp$iw),1,pk)
gha<-1.41*gforcedfree(vegp$lw*0.71,uz,tair,5,pk,5)
PAR<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=0.2)
gC<-layercond(PAR,vegp$gsmax,vegp$q50)
gv<-1/(1/gC+1/gha)
gtz<-phair(tair)*uz
gL<-1/(1/gha+1/gtz)
# Calculate absorbed radiation
Rabs<-.leafabs2(climdata$swrad,tme,tair,tground,lat,long,PAIt,PAIu,pLAI,vegp$x,vegp$refls,vegp$refw,
vegp$refg,vegp$vegem,climdata$skyem,dp,vegp$clump)
Th<-tleafS(tair,tground,relhum,pk,theta,gtt,gt0,gha,gv,gL,Rabs,vegp$vegem,soilp$b,
soilp$psi_e,soilp$Smax,surfwet,leafdens)
tz<-Th$tn
tleaf<-Th$tleaf
rh<-Th$rh
# Radiation
Rsw<-cansw(climdata$swrad,dp,tme=tme,lat=lat,long=long,x=vegp$x,l=PAIu,ref=vegp$refls)
# Calculate diffuse proportion
Rdif<-cantransdif(l=PAIu,ref=vegp$refls)*climdata$swrad*dp
dp<-Rdif/Rsw
dp[is.na(dp)]<-1
dp[is.infinite(dp)]<-1
dp[dp<0]<-0
dp[dp>1]<-1
Rlw<-canlw(tair,PAIu,1-vegp$vegem,climdata$skyem,vegp$clump)$lwin
}
metout<-data.frame(obs_time=climdata$obs_time,Tref=climdata$temp,Tloc=tz,tleaf=tleaf,
RHref=relhum,RHloc=rh,RSWloc=Rsw,RLWloc=Rlw,windspeed=uz,dp=dp)
# Consider snow
snow<-nmrout$snow
if (class(snow)[1] == "data.frame") {
snowtemp<-nmrout$SN1
} else snowtemp<-rep(0,length(tz))
mo<-nmrout$metout
snowdep<-mo$SNOWDEP
if (max(snowdep) > 0) {
mos <- .runmodelsnow(climdata,vegp,soilp,nmrout,reqhgt,lat,long,metopen,windhgt)
sel<-which(snowdep>0)
metout$Tloc[sel]<-mos$Tloc[sel]
metout$tleaf[sel]<-mos$tleaf[sel]
metout$RHloc[sel]<-mos$RHloc[sel]
metout$RSWloc[sel]<-mos$RSWloc[sel]
metout$RLWloc[sel]<-mos$RLWloc[sel]
metout$windspeed[sel]<-mos$windspeed[sel]
metout$dp[sel]<-mos$dp[sel]
}
return(metout)
}
#' Runs microclimate model in hourly timesteps with NicheMapR
#'
#' @description This function is a wrapper function for [runNMR()] and [runmodelS()] enabling
#' full microclimate model at any height above or below ground rapidly. NicheMapR is to
#' compute snowdepths, soil moistures and below-ground soil temperatures. The function
#' runmodelS is used to compute below or above canopy air temperatures.
#'
#' @param climdata data.frame of climate variables needed to run the run the model. The dataset should follow format
#' of [weather()]). Times must be in UTC.
#' @param prec vector of hourly or daily precipitation (mm) (see details).
#' @param vegp a list of vegetation parameters as returned by [microctools::habitatvars()].
#' @param soilp a list of soil parameters as returned by [soilinit()].
#' @param reqhgt height (m) for which microclimate is needed.
#' @param lat latitude of location (decimal degrees).
#' @param long longitude of location (decimal degrees).
#' @param altt elevation of location (m).
#' @param slope slope of ground surface (decimal degrees).
#' @param aspect aspect of ground surface (decimal degrees, 0 = north).
#' @param metopen optional logical indicating whether the wind measurement used as an input to
#' the model is from a nearby weather station located in open ground (TRUE) or above the canopy
#' for which temperatures are modelled (FALSE - see details)
#' @param windhgt height above ground of wind measurement. If `metopen` is FALSE, must be above
#' canopy.
#' @param surfwet proportion of leaf surface acting as free water surface.
#' @param groundem thermal emissivity of ground layer.
#' @param ERR Integrator error tolerance for soil temperature calculations.
#' @param cap Is organic cap present on soil surface? (1 = Yes, 0 = No).
#' @param hori Horizon angles (degrees), from 0 degrees azimuth (north) clockwise in 10 degree intervals.
#' @param maxpool Max depth for water pooling on the surface (mm), to account for runoff.
#' @param rainmult Rain multiplier for surface soil moisture (used to induce runon).
#' @param SoilMoist_Init Initial volumetric soil water content at each soil node (m^3/m^3)
#' @return a list of two objects:
#' @return a data.frame of microclimate variables for height `reqhgt` as returned by [runmodelS()]
#' @return a list of NicheMapR outputs as returned by [runNMR()]
#'
#' @import microctools
#' @import NicheMapR
#' @export
#' @seealso [runmodelS()], [runNMR()]
#'
#' @details This is a wrapper function for [runmodelS()] and [runNMR()] that allows fairly rapid
#' calaulation of microclimatic conditions below ground and/or above canopy. It uses the NicheMapR package
#' to calculate soil moisture-dependent soil temperatures, treating the vegetation canopy as a
#' single layer and then calculates air and leaf temperatures,air humidity and radiation using
#' runmodelS. Precipitation values are needed for computation of soil moisture and can either be
#' supplied as a vector of hourly values (must have identical number of values as in `climdata`) or
#' daily values (`climdata` must contain climate for entire days). It compares the length of the vector
#' to `climdata` to determine whether values are daily or hourly. In this function constant soil
#' poperties are are assumed throughout the soil profile. For more flexible implementation, [runNMR()]
#' and [runmodelS()] can be run seperately.
#'
#' @examples
#' require(NicheMapR)
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Run model using inbuilt weather datasets ~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' tme<-as.POSIXlt(weather$obs_time,tz="UTC")
#' vegp <- habitatvars(4, 50.2178, -5.32656, tme, m = 10) # Decidious broadleaf forest
#' soilp<- soilinit("Loam")
#' climrun <- runwithNMR(weather,dailyprecip,vegp,soilp,1,50.2178,-5.32656)
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Extract and view data from outputs ~~~~~~~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' metout<-climrun$metout
#' nmrout<-climrun$nmrout
#' soiltemps<-as.data.frame(nmrout$soil)
#' head(metout)
#' head(soiltemps)
#' head(nmrout$soilmoist)
#' head(nmrout$sunsnow)
#' head(nmrout$plant)
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Plot soil and air temperatures ~~~~~~~~~~~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' tsoil1<-soiltemps$D2.5cm
#' tsoil2<-soiltemps$D50cm
#' tair<-metout$Tloc
#' tref<-metout$Tref
#' tleaf<-metout$tleaf
#' Month<-as.POSIXct(metout$obs_time,tz="UTC")
#' tmn<-min(tsoil1,tsoil2,tair,tref,tleaf)
#' tmx<-max(tsoil1,tsoil2,tair,tref,tleaf)
#' par(mfrow=c(2,1))
#' # Air
#' plot(tref~Month,type="l",col=rgb(0,0,0,0.3),ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tleaf~Month,type="l",col=rgb(0,1,0,0.3),xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tair~Month,type="l",col=rgb(1,0,0,0.3),xlab="",ylab="",ylim=c(tmn,tmx))
#' # Soil
#' plot(tref~Month,type="l",col=rgb(0,0,0,0.5),ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil1~Month,type="l",col="red",xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil2~Month,type="l",col="blue",lwd=2,xlab="",ylab="",ylim=c(tmn,tmx))
#' # ====================================================================== #
#' # ~~~~~~~~~~~~~~~ Plot data by monthly mean ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
#' # ====================================================================== #
#' tme<-as.POSIXlt(metout$obs_time,tz="UTC")
#' tsoil1m<-aggregate(tsoil1,by=list(tme$mon),mean)$x
#' tsoil2m<-aggregate(tsoil2,by=list(tme$mon),mean)$x
#' tairm<-aggregate(tair,by=list(tme$mon),mean)$x
#' trefm<-aggregate(tref,by=list(tme$mon),mean)$x
#' tleafm<-aggregate(tleaf,by=list(tme$mon),mean)$x
#' tmn<-min(tsoil1m,tsoil2m,tairm,trefm,tleafm)
#' tmx<-max(tsoil1m,tsoil2m,tairm,trefm,tleafm)
#' Month<-c(1:12)
#' # Air
#' par(mfrow=c(2,1))
#' plot(trefm~Month,type="l",lwd=2,col="darkgray",ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tleafm~Month,type="l",lwd=2,col="darkgreen",xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tairm~Month,type="l",lwd=2,col="red",xlab="",ylab="",ylim=c(tmn,tmx))
# Soil
#' plot(trefm~Month,type="l",lwd=2,col="darkgray",ylab="Temperature",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil1m~Month,type="l",lwd=2,col="red",xlab="",ylab="",ylim=c(tmn,tmx))
#' par(new=T)
#' plot(tsoil2m~Month,type="l",lwd=2,col="blue",xlab="",ylab="",ylim=c(tmn,tmx))
#' # ====================================================================== #
runwithNMR <- function(climdata, prec, vegp, soilp, reqhgt, lat, long, altt = 0, slope = 0,
aspect = 0, metopen = TRUE, windhgt = 2, surfwet = 1, groundem = 0.95,
ERR = 1.5, cap = 1, hori = rep(0,36), maxpool =1000, rainmult = 1,
SoilMoist_Init = c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3),
animal = FALSE) {
if (!metopen & windhgt < vegp$hgt) {
stop("When metopen is FALSE, windhgt must be above canopy\n")
}
# (1) Unpack variables
tair<-climdata$temp
relhum<-climdata$relhum
pk<-climdata$pres
u<-climdata$windspeed
hgt<-vegp$hgt
# (1) Run NicheMapR
PAIt<-apply(vegp$PAI,2,sum)
LAI<-(vegp$pLAI*vegp$PAI)
LAI<-apply(LAI,2,mean)
pLAI<-LAI/PAIt
LREFL<-mean(pLAI*vegp$refls+(1-pLAI)*vegp$reflp)
soiltype<-as.character(soilp$Soil.type)
DEP<-c(0,2.5,5,10,15,20,30,50,100,200)
PE = rep(1.1, 19)
KS = rep(0.0037, 19)
BB = rep(4.5, 19)
BD = rep(1.3, 19)
DD = rep(2.65, 19)
nmrout<-runNMR(climdata,prec,lat,long,1.95,hgt,2,PAIt,vegp$x,pLAI,
vegp$clump,vegp$refg,LREFL,groundem,DEP,altt,slope,aspect,
ERR,soiltype,PE,KS,BB,BD,DD,cap,hori,maxpool,rainmult,SoilMoist_Init,
animal=animal)
if (reqhgt < 0) {
dep<-c(0,2.5,5,10,15,20,30,50,100,200)
soilz<- -reqhgt*100
dif<-abs(soilz-dep)
o<-order(dif)
dif1<-dif[o[1]]
dif2<-dif[o[2]]
tz1<-soilt[,o[1]+1]
tz2<-soilt[,o[2]+1]
tz<-tz1*(dif2/(dif1+dif2))+tz2*(dif1/(dif1+dif2))
metout <- data.frame(obs_time=climdata$obs_time,Tref=tair,Tloc=tz,
tleaf=-999,RHref=relhum,RHloc=-999)
} else if (reqhgt < (hgt+2)) {
metout <- runmodelS(climdata,vegp,soilp,nmrout,reqhgt,lat,long,metopen,windhgt,surfwet,groundem)
} else {
metout <- data.frame(obs_time=climdata$obs_time,Tref=tair,Tloc=tair,
tleaf=-999,RHref=relhum,RHloc=relhum)
warning("Height out of range. Output climate identical to input")
}
climdata$temp <- metout$Tloc
climdata$relhum <- metout$RHloc
climdata$swrad <- metout$RSWloc
climdata$windspeed <- metout$windspeed
nmrout2<-runNMR(climdata=climdata,prec=prec,lat=lat,long=long,
Usrhyt=1.95,Veghyt=hgt,Refhyt=reqhgt,PAIt,vegp$x,pLAI,
vegp$clump,vegp$refg,LREFL,groundem,DEP,altt,slope,aspect,
ERR,soiltype,PE,KS,BB,BD,DD,cap,hori,maxpool,rainmult,SoilMoist_Init,
animal=animal)
nmrouts<-nmrout2$nmrout
nmrouts$soil<-nmrout$nmrout$soil
nmrouts$soilmoist<-nmrout$nmrout$moist
return(list(metout = metout, nmrout = nmrouts))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.