#' @export
metGenRun <- function() {
mgcheckVariables()
## makeOutputNetCDF
makeNetcdfOut()
## DEFINE OUTPUT ARRAY
outData <- NULL
for (var in names(metGen$settings$outVars)) {
outData[[var]] <- array(NA, dim = c(metGen$settings$nx, metGen$settings$ny, metGen$derived$nOutStepDay))
}
nInStep <- metGen$derived$nInStepDay
nOutStep <- metGen$derived$nOutStepDay
maxStep <- max(nInStep, nOutStep)
if (nOutStep > nInStep) {
inrecs <- rep(1:nInStep, each = ceiling(nOutStep/nInStep))
outrecs <- rep(1:nOutStep)
} else {
inrecs <- rep(1:nInStep)
outrecs <- rep(1:nOutStep, each = ceiling(nInStep/nOutStep))
}
### THE MAIN LOOP
profile<-NULL
profile$start.time.total <- Sys.time()
for (iday in 1:metGen$derived$nday) {
metGen$current$timestep <- (metGen$derived$nInStepDay * (iday-1)) + 1
yday <- metGen$derived$inYDays[iday]
printf("Day: %d, date: %s\n", iday, metGen$derived$inDates[metGen$current$timestep])
## LOAD WHOLE DOMAIN FROM NETCDF
profile$start.time.read <- Sys.time()
inData <- readAllForcing(metGen$derived$inDates[metGen$current$timestep])
profile$end.time.read <- Sys.time()
# /*************************************************
# Precipitation
# *************************************************/
if(!is.null(outData$pr)) {
if(metGen$metadata$inVars$pr$enabled) { ## if precip is provided
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$pr[, ,outrecs[i]] <- inData$pr[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$pr[] <- 0
for(i in 1:maxStep) outData$pr[, , outrecs[i]] <- outData$pr[, , outrecs[i]] + ( inData$pr[, , inrecs[i]] / (nInStep/nOutStep) )
}
} else { ## if rainf and snowf are provided
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$pr[, ,outrecs[i]] <- inData$rainf[, , inrecs[i]] + inData$snowf[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$pr[] <- 0
for(i in 1:maxStep) outData$pr[, , outrecs[i]] <- outData$pr[, , outrecs[i]] + ( (inData$rainf[, , inrecs[i]] + inData$snowf[, , inrecs[i]]) / (nInStep/nOutStep) )
}
}
}
# /*************************************************
# Shortwave radiation
# *************************************************/
if (!is.null(outData$swdown)) {
if(metGen$metadata$inVars$swdown$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
radfrac <- rad_map_final_cr(metGen$derived$nOutStepDay, yday, gmt_float = 0, metGen$settings$lonlatbox)
swdown_day <- apply(inData$swdown, c(1,2), mean)
for(i in 1:maxStep) outData$swdown[, , outrecs[i]] <- radfrac[ , , outrecs[i]] * swdown_day
} else { ## aggregate to lower number of timesteps
outData$swdown[]<-0
for(i in 1:maxStep) outData$swdown[, , outrecs[i]] <- outData$swdown[, , outrecs[i]] + ( inData$swdown[, , inrecs[i]] / (nInStep/nOutStep) )
}
}
}
# /*************************************************
# Longwave radiation
# *************************************************/
if (!is.null(outData$lwdown)) {
if(metGen$metadata$inVars$lwdown$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$lwdown[, ,outrecs[i]] <- inData$lwdown[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$lwdown[] <- 0
for(i in 1:maxStep) outData$lwdown[, , outrecs[i]] <- outData$lwdown[, , outrecs[i]] + ( inData$lwdown[, , inrecs[i]] / (nInStep/nOutStep) )
}
}
}
# /*************************************************
# Wind
# *************************************************/
if (!is.null(outData$wind)) {
if(metGen$metadata$inVars$wind$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$wind[, ,outrecs[i]] <- inData$wind[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$wind[] <- 0
for(i in 1:maxStep) outData$wind[, , outrecs[i]] <- outData$wind[, , outrecs[i]] + ( inData$wind[, , inrecs[i]] / (nInStep/nOutStep) )
}
}
}
# /**************************************
# Atmospheric Pressure (Pa)
# **************************************/
if (!is.null(outData$psurf)) {
if(metGen$metadata$inVars$psurf$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$psurf[, ,outrecs[i]] <- inData$psurf[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$psurf[] <- 0
for(i in 1:maxStep) outData$psurf[, , outrecs[i]] <- outData$psurf[, , outrecs[i]] + ( inData$psurf[, , inrecs[i]] / (nInStep/nOutStep) )
}
}
}
# /**************************************
# Temperature
# **************************************/
if (!is.null(outData$tas)) {
if(metGen$metadata$inVars$tasmin$enabled && metGen$metadata$inVars$tasmax$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
if (nInStep > 1) stop(printf("Dissagregation of \"tasmin\" and \"tasmax\" into \"tas\" is only possible for daily input!"))
outData$tas <- set_max_min_lonlat_cr(inData$tasmin[,,1], inData$tasmax[,,1], yday, metGen$derived$nOutStepDay, metGen$settings$lonlatbox)
} else { ## aggregate to lower number of timesteps
outData$tas[] <- 0
for(i in 1:maxStep) outData$tas[, , outrecs[i]] <- outData$tas[, , outrecs[i]] +
( ( inData$tasmin[, , inrecs[i]] + inData$tasmax[, , inrecs[i]] ) / ((nInStep/nOutStep)*2) )
}
} else { ## if tas in provided as input
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$tas[, ,outrecs[i]] <- inData$tas[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$tas[] <- 0
for(i in 1:maxStep) outData$tas[, , outrecs[i]] <- outData$tas[, , outrecs[i]] + ( inData$tas[, , inrecs[i]] / (nInStep/nOutStep) )
}
}
}
# /**************************************
# Min and Max Temperature
# **************************************/
if (!is.null(outData$tasmin)) {
if(metGen$metadata$inVars$tas$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
stop(printf("disaggregation of \"tas\" into \"tasmin\" is not possible!"))
} else { ## aggregate to lower number of timesteps
outData$tasmin[] <- 0
for(i in 1:maxStep) outData$tasmin[, , outrecs[i]] <- pmin(outData$tasmin[, , outrecs[i]], inData$tas[, , inrecs[i]])
}
} else stop(printf("tas should be provided!\n"))
}
if (!is.null(outData$tasmax)) {
if(metGen$metadata$inVars$tas$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
stop(printf("disaggregation of \"tas\" into \"tasmax\" is not possible!"))
} else { ## aggregate to lower number of timesteps
outData$tasmax[] <- 0
for(i in 1:maxStep) outData$tasmax[, , outrecs[i]] <- pmax(outData$tasmax[, , outrecs[i]], inData$tas[, , inrecs[i]])
}
} else stop(printf("tas should be provided!\n"))
}
# /*************************************************
# Vapor pressure
# *************************************************/
# if(!is.null(metGen$settings$inVars$relhum) && !is.null(outData$vp)) {
# for(rec in 1:metGen$derived$nOutStepDay) {
# outData$vp <- set_vp_cr(outData$tas, inData$relhum[,,1], metGen$settings$nx, metGen$settings$ny, metGen$derived$nOutStepDay)
# }
# }
# if(!is.null(metGen$settings$inVars$qair) && !is.null(outData$vp)) {
# for(rec in 1:metGen$derived$nOutStepDay) {
# # outData$vp[,,rec] <- inData$qair[,,1] * inData$psurf[,,1] / metGen$constants$EPS
# # outData$vp[,,rec] <- mg_sh2vp(inData$qair[,,1], outData$tas[,,rec], inData$psurf[,,1])
# # outData$vp[,,rec] <- sh2vp(inData$qair[,,1], inData$psurf[,,1])
# outData$vp[,,rec] <- sh2vp(inData$qair[,,1], inData$psurf[,,1])
# }
# }
if (!is.null(outData$vp)) {
if (metGen$metadata$inVars$relhum$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) {
outData$vp[, , ] <- set_vp_cr(outData$tas[, , ], inData$relhum[, , inrecs[i]], metGen$settings$nx, metGen$settings$ny, metGen$derived$nOutStepDay)
}
} else { ## aggregate to lower number of timesteps
outData$vp[] <- 0
for(i in 1:maxStep) outData$vp[, , outrecs[i]] <- outData$vp[, , outrecs[i]] +
( set_vp_cr(outData$tas[, ,outrecs[i]], inData$relhum[, , inrecs[i]], metGen$settings$nx, metGen$settings$ny, metGen$derived$nOutStepDay) /
(nInStep/nOutStep) )
}
} else if(metGen$metadata$inVars$qair$enabled && metGen$metadata$inVars$psurf$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$vp[,,outrecs[i]] <- sh2vp(inData$qair[,,inrecs[i]], inData$psurf[,,inrecs[i]])
} else { ## aggregate to lower number of timesteps
outData$vp[] <- 0
for(i in 1:maxStep) outData$vp[, , outrecs[i]] <- outData$vp[, , outrecs[i]] + ( sh2vp(inData$qair[,,inrecs[i]], inData$psurf[,,inrecs[i]]) / (nInStep/nOutStep) )
}
} else if(metGen$metadata$inVars$vp$enabled) {
if (nInStep < nOutStep) { ## disaggregate to higher number of timesteps
for(i in 1:maxStep) outData$vp[, ,outrecs[i]] <- inData$vp[, , inrecs[i]]
} else { ## aggregate to lower number of timesteps
outData$vp[] <- 0
for(i in 1:maxStep) outData$vp[, , outrecs[i]] <- outData$vp[, , outrecs[i]] + ( inData$vp[, , inrecs[i]] / (nInStep/nOutStep) )
}
} else {
if (nInStep > 1) stop(printf("Cannot convert to vp!"))
}
}
## Convert to desired output unit
for (var in names(metGen$settings$outVars)) {
outData[[var]] <- convertUnit(outData[[var]],
metGen$metadata$outVars[[var]]$internal_units,
metGen$metadata$outVars[[var]]$output_units,
metGen$settings$outDt)
}
## ADD OUTPUT TO NETCDF
for (var in names(metGen$settings$outVars)) {
timeIndex <- metGen$derived$nOutStepDay*(iday-1)+1
metGen$settings$outVars[[var]]$ncid <- nc_open(metGen$settings$outVars[[var]]$filename, write = TRUE)
ncvar_put(metGen$settings$outVars[[var]]$ncid,
var,
outData[[var]][,,],
start = c(1, 1, timeIndex),
count = c(metGen$settings$nx, metGen$settings$ny, metGen$derived$nOutStepDay)
)
nc_close(metGen$settings$outVars[[var]]$ncid)
}
rm(inData)
}
ncellsTotal <- metGen$settings$nx*metGen$settings$ny
profile$end.time.total <- Sys.time()
totTime<-as.numeric(profile$end.time.total - profile$start.time.total, units = "secs")
cat(sprintf(" Total: %.1f seconds for %d day(s). So: 100 years will take: %.2f days\n",
totTime,
metGen$derived$nday,
(totTime * 365.25 * 100) / metGen$derived$nday / 86400
))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.