Nothing
R/thornthwaite.R
In ClimClass: Climate Classification According to Several Indices
NULL
#' @description Calculates Thornthwaite and Mather's water balance from monthly series of precipitation and temperature. Aimed at a classification of a site's climate according to its water balance features.
#'
#' @param series the monthly series of temperature and precipitation.
#' @param latitude latitude of the station in degrees.
#' @param clim_norm climatic normals.
#' @param first.yr first year of the period over which water balance is calculated. Default is \code{NULL} (calculations start with the first year of the series).
#' @param last.yr last year of the period over which water balance is calculated. Default is \code{NULL} (calculations stop with the last year of the series).
#' @param quant vector of quantiles for which water balance has to be assessed. Default is: min, 10th, 25th 50th, 75th, 90th, max.
#' @param snow.init initial water equivalent for snowpack (mm). Default is 20.
#' @param Tsnow maximum temperature (monthly mean) for precipitation to be treated as snowfall. Default is -1 degree C.
#' @param TAW maximum (field capacity) for soil water retention, and initial soil water content (mm). Default is 100.
#' @param fr.sn.acc fraction of snow that contributes to snowpack (0-1). 1 - fr.sn.acc is treated as liquid monthly precipitation Default is 0.95.
#' @param snow_melt_coeff monthly coefficient(s) for snowmelt. Default is 1.
#'
#'@import geosphere
#'
#'
#' @title Thornthwaite and Mather's water balance
#'
#' @author Giambattista Toller and Emanuele Eccel
#'
#' @return A \code{thornthwaite} S3 object, consisting on a list of two lists. The first (name: W_balance) is a list of data frames containing the monthly series of all indices, the second (name: quantiles) the relevant quantiles. See details for meanings of single variables.
#'
#' @details The algorithm for the calculation of water balance is adapted from Thornthwaite, 1948; Thornthwaite and Mather, 1955; Thornthwaite and Mather, 1957.
#'
#' \code{series} is a data frame with years, months, temperature and precipitation values. Names in series columns must include: year, month, Tn and Tx (minimum and maximum temperatures, respectively) or, as an alternative, Tm (mean temperatures), and P (mandatory).
#'
#' \code{clim_norm} is a monthly data frame of climate normals, with column names: "P", "Tn", "Tx", "Tm" (precipitation, minimum, maximum and mean temperature, respectively). It can be the output of function \code{\link{climate}}. If \code{clim_norm} is not NULL, any missing value in the monthly series is substituted by the corresponding climatic value in \code{clim_norm}.
#'
#' At any winter season, the maximum monthly snowpack height is attained in the last month before "spring" conditions (\code{Tm} >= \code{Tsnow}), even if a month with Tm < Tsnow may occur later.
#'
#' \code{snow_melt_coeff} is (are) the coefficient(s) for snow melt fraction(s) at any month where the condition for melting exists. If \code{snow_melt_coeff} = 1 (default), all the melting occurs in the first month when \code{Tm >= Tsnow}; if it is a vector, melting is spread over more than one month. If the sum of coefficients is less than 1, the residual melting occurs in one further month.
#'
#' The output function is a list of two lists of data frames (balance and quantile). In both lists, data frame (and names) are the following (all variables in mm):
#'
#' \code{Precipitation} (repeats input values);
#'
#' \code{Et0} (potential evapotranspiration);
#'
#' \code{Storage} (water stored in soil);
#'
#' \code{Prec. - Evap.} (difference between precipitation and potential evapotranspiration);
#'
#' \code{Deficit} (difference between potential and real evapotranspiration, due to water unavailability in soil);
#'
#' \code{Surplus} (water surplus in soil, routed to runoff).
#'
#' Please, refer to the quoted references for details.
#'
#' This function requires the function \code{\link{daylength}} (libr. \code{\link{geosphere}}).
#' @export
#'
#' @references
#' Thornthwaite, C. W., 1948: An Approach toward a Rational Classification of Climate. Geographical Review, Vol. 38, No. 1(Jan.):55-94.
#'
#' Thornthwaite, C. W., and Mather, J.R., 1955: The water balance. Publications in Climatology, Volume 8(1), Laboratory of Climatology
#'
#' Thornthwaite, C. W., and Mather, J.R., 1957: Instructions and tables for computing potential evapotranspiration and the water balance. Publications in climatology, Volume 10(3), Laboratory of Climatology
#'
#'
#' @examples
#'
#' data(Trent_climate)
#'
#'
#' # lista_cli is a list of data frames of the type "series",
#' # each one referring to one station - see function "climate".
#' # clima_81_10 is a list of data frames having climatic means
#' # of temperature and precipitation, each one referring to one station.
#' # It can be the output of function "climate".
#' library(geosphere) # required for function daylength
#' thornt_lst<-NULL
#' lista_cli <- lista_cli[1:3] ## lista_cli is reduced to diminish elapsed time of execution!
#' for(k in 1 : length(lista_cli[1:3])) {
#' thornt_lst[[k]]<-thornthwaite(series=lista_cli[[k]],
#' clim_norm=clima_81_10[[k]],
#' latitude = 46, first.yr=1981,
#' last.yr=2010, snow_melt_coeff=c(0.5,0.5 ) )
#' }
#' names(thornt_lst)<-names(lista_cli)
#'
#' # splits list into two lists
#' W_balance<-NULL; quantiles<-NULL
#' for(k in 1 : length(lista_cli))
#' {
#' W_balance[[k]]<-thornt_lst[[k]]$W_balance
#' quantiles[[k]]<-thornt_lst[[k]]$quantiles
#' }
#' names(W_balance)<-names(thornt_lst); names(quantiles)<-names(thornt_lst)
#'
#' @seealso \code{\link{climate}}, \code{\link{ExAtRa}}, \code{\link{plot.thornthwaite}}
thornthwaite<-function (series, latitude, clim_norm = NULL, first.yr = NULL,
last.yr = NULL, quant = c(0, 0.1, 0.25, 0.5, 0.75, 0.9, 1),
snow.init = 20, Tsnow = -1, TAW = 100, fr.sn.acc = 0.95,
snow_melt_coeff = 1)
{
P.exists<-"P" %in% names(series); Tn.exists<-"Tn" %in% names(series);Tx.exists<-"Tx" %in% names(series);Tm.exists<-"Tm" %in% names(series)
if(((!Tn.exists | !Tx.exists) & !Tm.exists) | !P.exists)
print("Both T and P series required!", quote=FALSE) else {
month_names <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
month_lengths <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31,
30, 31)
m <- c("min", paste(quant[-c(1, length(quant))] * 100, "%",
sep = ""), "max")
yDay <- strptime(paste(15, "/", 1:12, "/", 2014, sep = ""),
format = "%d/%m/%Y")$yday + 1
rel_daylength <- daylength(latitude, yDay)/12
chg <- ExAtRa(DOY = yDay, latitude = latitude, unit = "mm")
fTAW <- (-1.05e-07 * TAW^2 + 5.4225e-05 * TAW - 0.00878)
ultStorPrA <- TAW
if (is.null(first.yr))
first.yr <- min(series$year) else first.yr <- max(first.yr, min(series$year))
if (is.null(last.yr))
last.yr <- max(series$year) else last.yr <- min(last.yr, max(series$year))
years <- first.yr:last.yr
tPrec <- seq(1:12)
tET0 <- seq(1:12)
tStor <- seq(1:12)
tPmE <- seq(1:12)
tDef <- seq(1:12)
tSur <- seq(1:12)
snowpack_jan <- snow.init
for (ii in 1:length(years)) {
zzz <- data.frame(series[series$year == years[ii], ])
if(Tn.exists & Tx.exists & !Tm.exists)
zzz <- data.frame(zzz$P, (zzz$Tn + zzz$Tx)/2) else zzz <- data.frame(zzz$P, zzz$Tm)
names(zzz) <- c("Prec1", "Tmean1")
if (!is.null(clim_norm)) {
zzz$Prec1[is.na(zzz$Prec1)] <- clim_norm$P[is.na(zzz$Prec1)]
zzz$Tmean1[is.na(zzz$Tmean1)] <- clim_norm$Tm[is.na(zzz$Tmean1)]
}
if (sum(is.na(zzz$Prec1)) != 0 | sum(is.na(zzz$Tmean1)) != 0)
print(paste("Monthly values missing in year", first.yr +
ii - 1), quote = FALSE)
else {
Prec <- zzz$Prec1
Tmean <- zzz$Tmean1
Storage <- rep(TAW, length(Tmean))
Percola <- rep(0, length(Tmean))
TmeanCor <- Tmean
TmeanCor[TmeanCor < 0] <- 0
ITM <- ((TmeanCor/5)^1.514)
ITA <- sum(ITM)
exp <- 6.75e-07 * ITA^3 - 7.71e-05 * ITA^2 + 0.01792 *
ITA + 0.492
PET <- (16 * (10 * TmeanCor/ITA)^exp) * rel_daylength
SnowAcc <- rep(0, 12)
SnowPrec <- Prec
SnowPrec[Tmean >= Tsnow] <- 0
month_max <- max(min(which(Tmean >= Tsnow) - 1),
1)
SnowAcc_wint <- NULL
for (i in 1:month_max) {
if (i == 1)
SnowAcc_wint[i] <- snowpack_jan + SnowPrec[i] *
fr.sn.acc
else SnowAcc_wint[i] <- SnowAcc_wint[i - 1] +
SnowPrec[i] * fr.sn.acc
}
snowpack_ref <- SnowAcc_wint[month_max]
snow_depl <- NULL
SnowAcc[1] <- SnowAcc_wint[1]
snow_depl[1] <- SnowPrec[1] * (1 - fr.sn.acc)
if (Tmean[1] >= Tsnow) {
snow_depl[1] <- snow_melt_coeff[1] * SnowAcc[1]
SnowAcc[1] <- SnowAcc[1] - snow_depl[1]
}
count_thaw <- 0
for (i in 2:12) {
snow_depl[i] <- SnowPrec[i] * (1 - fr.sn.acc)
SnowAcc[i] <- SnowAcc[i - 1] + SnowPrec[i] *
fr.sn.acc
if (Tmean[i] >= Tsnow) {
count_thaw <- count_thaw + 1
if (count_thaw > length(snow_melt_coeff)) {
snow_depl[i] <- SnowAcc[i]
}
else snow_depl[i] <- snow_depl[i] + SnowAcc[i] *
snow_melt_coeff[count_thaw]
}
SnowAcc[i] <- SnowAcc[i] - snow_depl[i]
}
snowpack_jan <- SnowAcc[12]
Liq_Prec <- Prec
Liq_Prec[Tmean < Tsnow] <- Prec[Tmean < Tsnow] *
(1 - fr.sn.acc)
P.minus.ET <- Liq_Prec + snow_depl - PET
if (ii == 1)
last_Storage <- TAW
if (!is.na(P.minus.ET[1]) & P.minus.ET[1] > 0) {
Storage[1] <- last_Storage + P.minus.ET[1]
if (Storage[1] > TAW) {
Percola[1] <- Storage[1] - TAW
Storage[1] <- TAW
}
}
else {
PETvir <- (log10(TAW) - log10(last_Storage))/fTAW
Storage[1] <- TAW * 10^(-(PETvir + P.minus.ET[1]) *
fTAW)
}
for (i in 2:length(Storage)) {
if (!is.na(P.minus.ET[i]) & P.minus.ET[i] > 0) {
Storage[i] <- Storage[i - 1] + P.minus.ET[i]
if (Storage[i] > TAW) {
Percola[i] <- Storage[i] - TAW
Storage[i] <- TAW
}
}
else {
PETvir <- (log10(TAW) - log10(Storage[i - 1]))/fTAW
Storage[i] <- TAW * 10^(-(PETvir + P.minus.ET[i]) *
fTAW)
}
}
last_Storage <- Storage[12]
delta.sto <- c(0, diff(Storage))
ETr <- (Liq_Prec + snow_depl - delta.sto)
for (i in 1:length(ETr)) {
if (P.minus.ET[i] > 0)
ETr[i] <- PET[i]
}
Def <- PET - ETr
Def[Def < 0] <- 0
Sur <- Liq_Prec + snow_depl - ETr - (TAW - Storage)
Sur[Sur < 0] <- 0
tPrec <- data.frame(tPrec, Prec)
tET0 <- data.frame(tET0, PET)
tStor <- data.frame(tStor, Storage)
tPmE <- data.frame(tPmE, P.minus.ET)
tDef <- data.frame(tDef, Def)
tSur <- data.frame(tSur, Sur)
}
}
names(tPrec) <- c("Month", years)
row.names(tPrec) <- month_names
names(tET0) <- c("Month", years)
row.names(tET0) <- month_names
names(tStor) <- c("Month", years)
row.names(tStor) <- month_names
names(tPmE) <- c("Month", years)
row.names(tPmE) <- month_names
names(tDef) <- c("Month", years)
row.names(tDef) <- month_names
names(tSur) <- c("Month", years)
row.names(tSur) <- month_names
list_thornt <- list(round(tPrec[-1], 1), round(tET0[-1],
1), round(tStor[-1], 1), round(tPmE[-1], 1), round(tDef[-1],
1), round(tSur[-1], 1))
names(list_thornt) <- c("Precipitation", "Et0", "Storage",
"Prec. - Evap.", "Deficit", "Surplus")
qPrec <- data.frame(quant)
qET0 <- data.frame(quant)
qStor <- data.frame(quant)
qPmE <- data.frame(quant)
qDef <- data.frame(quant)
qSur <- data.frame(quant)
ttPrec <- t(tPrec)
ttET0 <- t(tET0)
ttStor <- t(tStor)
ttPmE <- t(tPmE)
ttDef <- t(tDef)
ttSur <- t(tSur)
qPrec <- as.data.frame(apply(ttPrec[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qPrec) <- month_names
qET0 <- as.data.frame(apply(ttET0[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qET0) <- month_names
qStor <- as.data.frame(apply(ttStor[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qStor) <- month_names
qPmE <- as.data.frame(apply(ttPmE[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qPmE) <- month_names
qDef <- as.data.frame(apply(ttDef[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qDef) <- month_names
qSur <- as.data.frame(apply(ttSur[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qSur) <- month_names
list_quant <- list(round(qPrec, 1), round(qET0, 1), round(qStor,
1), round(qPmE, 1), round(qDef, 1), round(qSur, 1))
names(list_quant) <- names(list_thornt)
list_2.lists <- list(W_balance = list_thornt, quantiles = list_quant)
class(list_2.lists) <- "thornthwaite"
return(list_2.lists)
} # END else (temperature and precipitation data do exist)
}
Try the ClimClass package in your browser
Any scripts or data that you put into this service are public.
ClimClass documentation built on May 2, 2019, 2:20 p.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.
NULL
#' @description Calculates Thornthwaite and Mather's water balance from monthly series of precipitation and temperature. Aimed at a classification of a site's climate according to its water balance features.
#'
#' @param series the monthly series of temperature and precipitation.
#' @param latitude latitude of the station in degrees.
#' @param clim_norm climatic normals.
#' @param first.yr first year of the period over which water balance is calculated. Default is \code{NULL} (calculations start with the first year of the series).
#' @param last.yr last year of the period over which water balance is calculated. Default is \code{NULL} (calculations stop with the last year of the series).
#' @param quant vector of quantiles for which water balance has to be assessed. Default is: min, 10th, 25th 50th, 75th, 90th, max.
#' @param snow.init initial water equivalent for snowpack (mm). Default is 20.
#' @param Tsnow maximum temperature (monthly mean) for precipitation to be treated as snowfall. Default is -1 degree C.
#' @param TAW maximum (field capacity) for soil water retention, and initial soil water content (mm). Default is 100.
#' @param fr.sn.acc fraction of snow that contributes to snowpack (0-1). 1 - fr.sn.acc is treated as liquid monthly precipitation Default is 0.95.
#' @param snow_melt_coeff monthly coefficient(s) for snowmelt. Default is 1.
#'
#'@import geosphere
#'
#'
#' @title Thornthwaite and Mather's water balance
#'
#' @author Giambattista Toller and Emanuele Eccel
#'
#' @return A \code{thornthwaite} S3 object, consisting on a list of two lists. The first (name: W_balance) is a list of data frames containing the monthly series of all indices, the second (name: quantiles) the relevant quantiles. See details for meanings of single variables.
#'
#' @details The algorithm for the calculation of water balance is adapted from Thornthwaite, 1948; Thornthwaite and Mather, 1955; Thornthwaite and Mather, 1957.
#'
#' \code{series} is a data frame with years, months, temperature and precipitation values. Names in series columns must include: year, month, Tn and Tx (minimum and maximum temperatures, respectively) or, as an alternative, Tm (mean temperatures), and P (mandatory).
#'
#' \code{clim_norm} is a monthly data frame of climate normals, with column names: "P", "Tn", "Tx", "Tm" (precipitation, minimum, maximum and mean temperature, respectively). It can be the output of function \code{\link{climate}}. If \code{clim_norm} is not NULL, any missing value in the monthly series is substituted by the corresponding climatic value in \code{clim_norm}.
#'
#' At any winter season, the maximum monthly snowpack height is attained in the last month before "spring" conditions (\code{Tm} >= \code{Tsnow}), even if a month with Tm < Tsnow may occur later.
#'
#' \code{snow_melt_coeff} is (are) the coefficient(s) for snow melt fraction(s) at any month where the condition for melting exists. If \code{snow_melt_coeff} = 1 (default), all the melting occurs in the first month when \code{Tm >= Tsnow}; if it is a vector, melting is spread over more than one month. If the sum of coefficients is less than 1, the residual melting occurs in one further month.
#'
#' The output function is a list of two lists of data frames (balance and quantile). In both lists, data frame (and names) are the following (all variables in mm):
#'
#' \code{Precipitation} (repeats input values);
#'
#' \code{Et0} (potential evapotranspiration);
#'
#' \code{Storage} (water stored in soil);
#'
#' \code{Prec. - Evap.} (difference between precipitation and potential evapotranspiration);
#'
#' \code{Deficit} (difference between potential and real evapotranspiration, due to water unavailability in soil);
#'
#' \code{Surplus} (water surplus in soil, routed to runoff).
#'
#' Please, refer to the quoted references for details.
#'
#' This function requires the function \code{\link{daylength}} (libr. \code{\link{geosphere}}).
#' @export
#'
#' @references
#' Thornthwaite, C. W., 1948: An Approach toward a Rational Classification of Climate. Geographical Review, Vol. 38, No. 1(Jan.):55-94.
#'
#' Thornthwaite, C. W., and Mather, J.R., 1955: The water balance. Publications in Climatology, Volume 8(1), Laboratory of Climatology
#'
#' Thornthwaite, C. W., and Mather, J.R., 1957: Instructions and tables for computing potential evapotranspiration and the water balance. Publications in climatology, Volume 10(3), Laboratory of Climatology
#'
#'
#' @examples
#'
#' data(Trent_climate)
#'
#'
#' # lista_cli is a list of data frames of the type "series",
#' # each one referring to one station - see function "climate".
#' # clima_81_10 is a list of data frames having climatic means
#' # of temperature and precipitation, each one referring to one station.
#' # It can be the output of function "climate".
#' library(geosphere) # required for function daylength
#' thornt_lst<-NULL
#' lista_cli <- lista_cli[1:3] ## lista_cli is reduced to diminish elapsed time of execution!
#' for(k in 1 : length(lista_cli[1:3])) {
#' thornt_lst[[k]]<-thornthwaite(series=lista_cli[[k]],
#' clim_norm=clima_81_10[[k]],
#' latitude = 46, first.yr=1981,
#' last.yr=2010, snow_melt_coeff=c(0.5,0.5 ) )
#' }
#' names(thornt_lst)<-names(lista_cli)
#'
#' # splits list into two lists
#' W_balance<-NULL; quantiles<-NULL
#' for(k in 1 : length(lista_cli))
#' {
#' W_balance[[k]]<-thornt_lst[[k]]$W_balance
#' quantiles[[k]]<-thornt_lst[[k]]$quantiles
#' }
#' names(W_balance)<-names(thornt_lst); names(quantiles)<-names(thornt_lst)
#'
#' @seealso \code{\link{climate}}, \code{\link{ExAtRa}}, \code{\link{plot.thornthwaite}}
thornthwaite<-function (series, latitude, clim_norm = NULL, first.yr = NULL,
last.yr = NULL, quant = c(0, 0.1, 0.25, 0.5, 0.75, 0.9, 1),
snow.init = 20, Tsnow = -1, TAW = 100, fr.sn.acc = 0.95,
snow_melt_coeff = 1)
{
P.exists<-"P" %in% names(series); Tn.exists<-"Tn" %in% names(series);Tx.exists<-"Tx" %in% names(series);Tm.exists<-"Tm" %in% names(series)
if(((!Tn.exists | !Tx.exists) & !Tm.exists) | !P.exists)
print("Both T and P series required!", quote=FALSE) else {
month_names <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
month_lengths <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31,
30, 31)
m <- c("min", paste(quant[-c(1, length(quant))] * 100, "%",
sep = ""), "max")
yDay <- strptime(paste(15, "/", 1:12, "/", 2014, sep = ""),
format = "%d/%m/%Y")$yday + 1
rel_daylength <- daylength(latitude, yDay)/12
chg <- ExAtRa(DOY = yDay, latitude = latitude, unit = "mm")
fTAW <- (-1.05e-07 * TAW^2 + 5.4225e-05 * TAW - 0.00878)
ultStorPrA <- TAW
if (is.null(first.yr))
first.yr <- min(series$year) else first.yr <- max(first.yr, min(series$year))
if (is.null(last.yr))
last.yr <- max(series$year) else last.yr <- min(last.yr, max(series$year))
years <- first.yr:last.yr
tPrec <- seq(1:12)
tET0 <- seq(1:12)
tStor <- seq(1:12)
tPmE <- seq(1:12)
tDef <- seq(1:12)
tSur <- seq(1:12)
snowpack_jan <- snow.init
for (ii in 1:length(years)) {
zzz <- data.frame(series[series$year == years[ii], ])
if(Tn.exists & Tx.exists & !Tm.exists)
zzz <- data.frame(zzz$P, (zzz$Tn + zzz$Tx)/2) else zzz <- data.frame(zzz$P, zzz$Tm)
names(zzz) <- c("Prec1", "Tmean1")
if (!is.null(clim_norm)) {
zzz$Prec1[is.na(zzz$Prec1)] <- clim_norm$P[is.na(zzz$Prec1)]
zzz$Tmean1[is.na(zzz$Tmean1)] <- clim_norm$Tm[is.na(zzz$Tmean1)]
}
if (sum(is.na(zzz$Prec1)) != 0 | sum(is.na(zzz$Tmean1)) != 0)
print(paste("Monthly values missing in year", first.yr +
ii - 1), quote = FALSE)
else {
Prec <- zzz$Prec1
Tmean <- zzz$Tmean1
Storage <- rep(TAW, length(Tmean))
Percola <- rep(0, length(Tmean))
TmeanCor <- Tmean
TmeanCor[TmeanCor < 0] <- 0
ITM <- ((TmeanCor/5)^1.514)
ITA <- sum(ITM)
exp <- 6.75e-07 * ITA^3 - 7.71e-05 * ITA^2 + 0.01792 *
ITA + 0.492
PET <- (16 * (10 * TmeanCor/ITA)^exp) * rel_daylength
SnowAcc <- rep(0, 12)
SnowPrec <- Prec
SnowPrec[Tmean >= Tsnow] <- 0
month_max <- max(min(which(Tmean >= Tsnow) - 1),
1)
SnowAcc_wint <- NULL
for (i in 1:month_max) {
if (i == 1)
SnowAcc_wint[i] <- snowpack_jan + SnowPrec[i] *
fr.sn.acc
else SnowAcc_wint[i] <- SnowAcc_wint[i - 1] +
SnowPrec[i] * fr.sn.acc
}
snowpack_ref <- SnowAcc_wint[month_max]
snow_depl <- NULL
SnowAcc[1] <- SnowAcc_wint[1]
snow_depl[1] <- SnowPrec[1] * (1 - fr.sn.acc)
if (Tmean[1] >= Tsnow) {
snow_depl[1] <- snow_melt_coeff[1] * SnowAcc[1]
SnowAcc[1] <- SnowAcc[1] - snow_depl[1]
}
count_thaw <- 0
for (i in 2:12) {
snow_depl[i] <- SnowPrec[i] * (1 - fr.sn.acc)
SnowAcc[i] <- SnowAcc[i - 1] + SnowPrec[i] *
fr.sn.acc
if (Tmean[i] >= Tsnow) {
count_thaw <- count_thaw + 1
if (count_thaw > length(snow_melt_coeff)) {
snow_depl[i] <- SnowAcc[i]
}
else snow_depl[i] <- snow_depl[i] + SnowAcc[i] *
snow_melt_coeff[count_thaw]
}
SnowAcc[i] <- SnowAcc[i] - snow_depl[i]
}
snowpack_jan <- SnowAcc[12]
Liq_Prec <- Prec
Liq_Prec[Tmean < Tsnow] <- Prec[Tmean < Tsnow] *
(1 - fr.sn.acc)
P.minus.ET <- Liq_Prec + snow_depl - PET
if (ii == 1)
last_Storage <- TAW
if (!is.na(P.minus.ET[1]) & P.minus.ET[1] > 0) {
Storage[1] <- last_Storage + P.minus.ET[1]
if (Storage[1] > TAW) {
Percola[1] <- Storage[1] - TAW
Storage[1] <- TAW
}
}
else {
PETvir <- (log10(TAW) - log10(last_Storage))/fTAW
Storage[1] <- TAW * 10^(-(PETvir + P.minus.ET[1]) *
fTAW)
}
for (i in 2:length(Storage)) {
if (!is.na(P.minus.ET[i]) & P.minus.ET[i] > 0) {
Storage[i] <- Storage[i - 1] + P.minus.ET[i]
if (Storage[i] > TAW) {
Percola[i] <- Storage[i] - TAW
Storage[i] <- TAW
}
}
else {
PETvir <- (log10(TAW) - log10(Storage[i - 1]))/fTAW
Storage[i] <- TAW * 10^(-(PETvir + P.minus.ET[i]) *
fTAW)
}
}
last_Storage <- Storage[12]
delta.sto <- c(0, diff(Storage))
ETr <- (Liq_Prec + snow_depl - delta.sto)
for (i in 1:length(ETr)) {
if (P.minus.ET[i] > 0)
ETr[i] <- PET[i]
}
Def <- PET - ETr
Def[Def < 0] <- 0
Sur <- Liq_Prec + snow_depl - ETr - (TAW - Storage)
Sur[Sur < 0] <- 0
tPrec <- data.frame(tPrec, Prec)
tET0 <- data.frame(tET0, PET)
tStor <- data.frame(tStor, Storage)
tPmE <- data.frame(tPmE, P.minus.ET)
tDef <- data.frame(tDef, Def)
tSur <- data.frame(tSur, Sur)
}
}
names(tPrec) <- c("Month", years)
row.names(tPrec) <- month_names
names(tET0) <- c("Month", years)
row.names(tET0) <- month_names
names(tStor) <- c("Month", years)
row.names(tStor) <- month_names
names(tPmE) <- c("Month", years)
row.names(tPmE) <- month_names
names(tDef) <- c("Month", years)
row.names(tDef) <- month_names
names(tSur) <- c("Month", years)
row.names(tSur) <- month_names
list_thornt <- list(round(tPrec[-1], 1), round(tET0[-1],
1), round(tStor[-1], 1), round(tPmE[-1], 1), round(tDef[-1],
1), round(tSur[-1], 1))
names(list_thornt) <- c("Precipitation", "Et0", "Storage",
"Prec. - Evap.", "Deficit", "Surplus")
qPrec <- data.frame(quant)
qET0 <- data.frame(quant)
qStor <- data.frame(quant)
qPmE <- data.frame(quant)
qDef <- data.frame(quant)
qSur <- data.frame(quant)
ttPrec <- t(tPrec)
ttET0 <- t(tET0)
ttStor <- t(tStor)
ttPmE <- t(tPmE)
ttDef <- t(tDef)
ttSur <- t(tSur)
qPrec <- as.data.frame(apply(ttPrec[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qPrec) <- month_names
qET0 <- as.data.frame(apply(ttET0[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qET0) <- month_names
qStor <- as.data.frame(apply(ttStor[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qStor) <- month_names
qPmE <- as.data.frame(apply(ttPmE[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qPmE) <- month_names
qDef <- as.data.frame(apply(ttDef[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qDef) <- month_names
qSur <- as.data.frame(apply(ttSur[-1, ], FUN = quantile,
probs = quant, na.rm = T, MARGIN = 2))
names(qSur) <- month_names
list_quant <- list(round(qPrec, 1), round(qET0, 1), round(qStor,
1), round(qPmE, 1), round(qDef, 1), round(qSur, 1))
names(list_quant) <- names(list_thornt)
list_2.lists <- list(W_balance = list_thornt, quantiles = list_quant)
class(list_2.lists) <- "thornthwaite"
return(list_2.lists)
} # END else (temperature and precipitation data do exist)
}
Try the ClimClass package in your browser
Any scripts or data that you put into this service are public.
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.