R/thornthwaite.R
In ecor/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.10,0.25,0.50,0.75,0.90,1.00),
snow.init=20, Tsnow=-1, TAW=100, fr.sn.acc=0.95, snow_melt_coeff=1)
{
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")
# creates a set of DOYs centered in each month
yDay <- strptime(paste(15,"/",1:12,"/",2014,sep=""), format="%d/%m/%Y")$yday + 1
rel_daylength<-daylength(latitude,yDay)/12
# calculates extra-atmospheric radiation for each day in the set
chg<- ExAtRa(DOY=yDay,latitude=latitude, unit="mm")
# calculates total available water
fTAW<-(-0.000000105*TAW^2+0.000054225*TAW-0.00878)
# first storage value = TAW
ultStorPrA<-TAW
# determines calculation period
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)
# creates a data frame
snowpack_jan<-snow.init
for(ii in 1:length(years))
{
zzz<-data.frame(series[series$year==years[ii],])
zzz<-data.frame(zzz[3], (zzz$Tn + zzz$Tx)/2, zzz$Tn, zzz$Tx); names(zzz)<-c("Prec1", "Tmean1","Tn1", "Tx1")
# substitutes missing values with those of clim_norm, if passed to function
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)]
zzz$Tn1[is.na(zzz$Tn1)]<-clim_norm$Tn[is.na(zzz$Tn1)]
zzz$Tx1[is.na(zzz$Tx1)]<-clim_norm$Tx[is.na(zzz$Tx1)]
}
# check on residual missing values
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[,1]
Tmean <-zzz[,2]
Tn <-zzz[,3]
Tx <-zzz[,4]
Storage<-rep(TAW,length(Tmean))
Percola<-rep(0,length(Tmean))
# calculation of thermal index
TmeanCor<-Tmean
TmeanCor[TmeanCor<0]<-0
ITM<-((TmeanCor/5)^1.514) # monthly
ITA<-sum(ITM) # annual
# calculation of PET
exp<- 0.000000675*ITA^3 - 0.0000771 *ITA^2 + 0.01792*ITA + 0.492
PET<-(16*(10*TmeanCor/ITA)^exp)*rel_daylength
# snow accumulation
SnowAcc<-rep(0,12)
SnowPrec<-Prec; SnowPrec[Tmean>=Tsnow]<-0
# calculates maximum snowpack before first thaw month
month_max<-max(min(which(Tmean>=Tsnow)-1), 1) # month before first with Tmean >= Tsnow
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) # 0 if Tmean>=Tsnow
if(Tmean[1]>=Tsnow)
{
snow_depl[1] <- snow_melt_coeff[1]*SnowAcc[1]
SnowAcc[1]<- SnowAcc[1]-snow_depl[1]
}
# following months
count_thaw<-0
for(i in 2:12)
{
snow_depl[i]<- SnowPrec[i]*(1-fr.sn.acc) # 0 if Tmean>=Tsnow
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]
# calculation of difference between Prec and PET
Liq_Prec<- Prec; Liq_Prec[Tmean<Tsnow] <- Prec[Tmean<Tsnow]*(1-fr.sn.acc)
P.minus.ET<- Liq_Prec + snow_depl - PET
# water balance
if(ii==1) # first year, field capacity
last_Storage<-TAW
# first month
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 # virtual PET
Storage[1]<-TAW*10^(-(PETvir + P.minus.ET[1])*fTAW)
}
# following months
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)
}
}
# saves the last storage
last_Storage<-Storage[12]
# calculates real ET as difference between prec. and storage
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]
}
# deficit: difference between real and potential ET (non-negative only)
Def <- PET - ETr; Def[Def<0]<-0
# water surplus
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)
} # end else (no missing monthly data)
} # end cycle on years
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")
# quantile tables
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)
}
ecor/ClimClass documentation built on May 15, 2019, 8 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.10,0.25,0.50,0.75,0.90,1.00),
snow.init=20, Tsnow=-1, TAW=100, fr.sn.acc=0.95, snow_melt_coeff=1)
{
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")
# creates a set of DOYs centered in each month
yDay <- strptime(paste(15,"/",1:12,"/",2014,sep=""), format="%d/%m/%Y")$yday + 1
rel_daylength<-daylength(latitude,yDay)/12
# calculates extra-atmospheric radiation for each day in the set
chg<- ExAtRa(DOY=yDay,latitude=latitude, unit="mm")
# calculates total available water
fTAW<-(-0.000000105*TAW^2+0.000054225*TAW-0.00878)
# first storage value = TAW
ultStorPrA<-TAW
# determines calculation period
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)
# creates a data frame
snowpack_jan<-snow.init
for(ii in 1:length(years))
{
zzz<-data.frame(series[series$year==years[ii],])
zzz<-data.frame(zzz[3], (zzz$Tn + zzz$Tx)/2, zzz$Tn, zzz$Tx); names(zzz)<-c("Prec1", "Tmean1","Tn1", "Tx1")
# substitutes missing values with those of clim_norm, if passed to function
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)]
zzz$Tn1[is.na(zzz$Tn1)]<-clim_norm$Tn[is.na(zzz$Tn1)]
zzz$Tx1[is.na(zzz$Tx1)]<-clim_norm$Tx[is.na(zzz$Tx1)]
}
# check on residual missing values
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[,1]
Tmean <-zzz[,2]
Tn <-zzz[,3]
Tx <-zzz[,4]
Storage<-rep(TAW,length(Tmean))
Percola<-rep(0,length(Tmean))
# calculation of thermal index
TmeanCor<-Tmean
TmeanCor[TmeanCor<0]<-0
ITM<-((TmeanCor/5)^1.514) # monthly
ITA<-sum(ITM) # annual
# calculation of PET
exp<- 0.000000675*ITA^3 - 0.0000771 *ITA^2 + 0.01792*ITA + 0.492
PET<-(16*(10*TmeanCor/ITA)^exp)*rel_daylength
# snow accumulation
SnowAcc<-rep(0,12)
SnowPrec<-Prec; SnowPrec[Tmean>=Tsnow]<-0
# calculates maximum snowpack before first thaw month
month_max<-max(min(which(Tmean>=Tsnow)-1), 1) # month before first with Tmean >= Tsnow
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) # 0 if Tmean>=Tsnow
if(Tmean[1]>=Tsnow)
{
snow_depl[1] <- snow_melt_coeff[1]*SnowAcc[1]
SnowAcc[1]<- SnowAcc[1]-snow_depl[1]
}
# following months
count_thaw<-0
for(i in 2:12)
{
snow_depl[i]<- SnowPrec[i]*(1-fr.sn.acc) # 0 if Tmean>=Tsnow
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]
# calculation of difference between Prec and PET
Liq_Prec<- Prec; Liq_Prec[Tmean<Tsnow] <- Prec[Tmean<Tsnow]*(1-fr.sn.acc)
P.minus.ET<- Liq_Prec + snow_depl - PET
# water balance
if(ii==1) # first year, field capacity
last_Storage<-TAW
# first month
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 # virtual PET
Storage[1]<-TAW*10^(-(PETvir + P.minus.ET[1])*fTAW)
}
# following months
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)
}
}
# saves the last storage
last_Storage<-Storage[12]
# calculates real ET as difference between prec. and storage
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]
}
# deficit: difference between real and potential ET (non-negative only)
Def <- PET - ETr; Def[Def<0]<-0
# water surplus
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)
} # end else (no missing monthly data)
} # end cycle on years
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")
# quantile tables
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)
}
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.