diff --git a/R/utils.R b/R/utils.R
index 4e961f9144a77095700d7d9b9dd508b459ae8122..a36aa1b1e22737806998118047014c3e9bc6ab33 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -177,7 +177,7 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
     }
   )
   names(lQsim) <- griwrm$id[QsimRows]
-  dfQsim <- cbind(data.frame(DatesR = as.POSIXct(InputsModel[[1]]$DatesR[IndPeriod_Run])),
+  dfQsim <- cbind(data.frame(DatesR = InputsModel[[1]]$DatesR[IndPeriod_Run]),
                   do.call(cbind,lQsim) / attr(InputsModel, "TimeStep"))
   class(dfQsim) <- c("Qm3s", class(dfQsim)) # For S3 methods
   return(dfQsim)
diff --git a/data-raw/Severn.R b/data-raw/Severn.R
index 5535ff91fc05663950a9fc17249e574eb632bdab..db720be1031c9481c2b04cb91611cba091d76c42 100644
--- a/data-raw/Severn.R
+++ b/data-raw/Severn.R
@@ -21,7 +21,7 @@ read_CAMELS_ts <- function(file, cut) {
   df <- read.csv(file)
   iCut <- which(df$date == cut)
   df <- df[iCut:nrow(df),]
-  df$DatesR <- as.POSIXct(df$date)
+  df$DatesR <- as.POSIXct(df$date, tz = "UTC")
 
   return(df[,c("DatesR", "precipitation", "peti", "discharge_spec")])
 }
diff --git a/data/Severn.rda b/data/Severn.rda
index 1d26cfd47d14b9716071c14f9346a4c05d61fbf5..1a6827656f3f7eeb2c655f4e5cc3cd5fdcab0aeb 100644
Binary files a/data/Severn.rda and b/data/Severn.rda differ
diff --git a/man-examples/RunModel.GRiwrmInputsModel.R b/man-examples/RunModel.GRiwrmInputsModel.R
index 9aafca1133e02849179b920c0f46fe4c3548c5c0..c84a96a8048755b34dc6a9e936fef0184114e70b 100644
--- a/man-examples/RunModel.GRiwrmInputsModel.R
+++ b/man-examples/RunModel.GRiwrmInputsModel.R
@@ -192,8 +192,8 @@ Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
 
 # Plot the diverted flow for the year 2003
 Ind_Plot <- which(
-  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01") &
-  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31")
+  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
diff --git a/man-examples/RunModel.Supervisor.R b/man-examples/RunModel.Supervisor.R
index d1aa0f1871052e07f55634bad8d67f126388d306..91bcf4fc1a46e1e509653478aae6c078a4d56559 100644
--- a/man-examples/RunModel.Supervisor.R
+++ b/man-examples/RunModel.Supervisor.R
@@ -45,7 +45,7 @@ PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
 Qobs <- data.frame(
   Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
 )
-Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01")), 1]
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
 
 # InputsModel object
 IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
@@ -96,7 +96,7 @@ CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN =
 
 # GRiwrmRunOptions object simulation of the hydrological year 2002-2003
 IndPeriod_Run <- seq.int(
-  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15")),
+  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15", tz = "UTC")),
   length.out = 365
 )
 IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
diff --git a/man/CreateController.Rd b/man/CreateController.Rd
index d474e737805e487df315f47ec164880a254e8cbd..470638b9a398cc89b6b837b0953dcac4c128b01d 100644
--- a/man/CreateController.Rd
+++ b/man/CreateController.Rd
@@ -90,7 +90,7 @@ PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
 Qobs <- data.frame(
   Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
 )
-Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01")), 1]
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
 
 # InputsModel object
 IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
@@ -141,7 +141,7 @@ CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN =
 
 # GRiwrmRunOptions object simulation of the hydrological year 2002-2003
 IndPeriod_Run <- seq.int(
-  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15")),
+  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15", tz = "UTC")),
   length.out = 365
 )
 IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
diff --git a/man/CreateInputsModel.GRiwrm.Rd b/man/CreateInputsModel.GRiwrm.Rd
index 608f953b2828bddf235714cd12848cb27a501696..9a438b0061505b4c1de25efb53f0d3191ffb97bc 100644
--- a/man/CreateInputsModel.GRiwrm.Rd
+++ b/man/CreateInputsModel.GRiwrm.Rd
@@ -295,8 +295,8 @@ Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
 
 # Plot the diverted flow for the year 2003
 Ind_Plot <- which(
-  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01") &
-  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31")
+  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
diff --git a/man/CreateRunOptions.Rd b/man/CreateRunOptions.Rd
index 56ed9eb1c341e5f00ed43cdb418f9b8b12b94510..51c0478531a6d9864f02d6f654f2b9dbce81ecdd 100644
--- a/man/CreateRunOptions.Rd
+++ b/man/CreateRunOptions.Rd
@@ -238,8 +238,8 @@ Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
 
 # Plot the diverted flow for the year 2003
 Ind_Plot <- which(
-  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01") &
-  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31")
+  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
diff --git a/man/CreateSupervisor.Rd b/man/CreateSupervisor.Rd
index d68a3c417056a7af209568e5a14cff412a8de83e..d7585184a226e42391be70e9d83b77b74a1b574a 100644
--- a/man/CreateSupervisor.Rd
+++ b/man/CreateSupervisor.Rd
@@ -72,7 +72,7 @@ PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
 Qobs <- data.frame(
   Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
 )
-Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01")), 1]
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
 
 # InputsModel object
 IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
@@ -123,7 +123,7 @@ CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN =
 
 # GRiwrmRunOptions object simulation of the hydrological year 2002-2003
 IndPeriod_Run <- seq.int(
-  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15")),
+  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15", tz = "UTC")),
   length.out = 365
 )
 IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
diff --git a/man/RunModel.GRiwrmInputsModel.Rd b/man/RunModel.GRiwrmInputsModel.Rd
index 097b74dc95bb08e264716e443825d597c1c22039..a87734567e37ae701cdfbea76d0696324bba1bf0 100644
--- a/man/RunModel.GRiwrmInputsModel.Rd
+++ b/man/RunModel.GRiwrmInputsModel.Rd
@@ -216,8 +216,8 @@ Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
 
 # Plot the diverted flow for the year 2003
 Ind_Plot <- which(
-  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01") &
-  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31")
+  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
diff --git a/man/RunModel.Supervisor.Rd b/man/RunModel.Supervisor.Rd
index 667be3f8e295c1037529dd946745d93b120c2261..3aa4f6b7f1045ffd363f5c09d4a80249192e246a 100644
--- a/man/RunModel.Supervisor.Rd
+++ b/man/RunModel.Supervisor.Rd
@@ -69,7 +69,7 @@ PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
 Qobs <- data.frame(
   Dam = BasinsObs$`54095`$discharge_spec * griwrm$area[griwrm$id == "54095"] * 1E3
 )
-Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01")), 1]
+Qobs[,] <- Qobs[which(DatesR == as.POSIXct("2002-10-01", tz = "UTC")), 1]
 
 # InputsModel object
 IM_severn <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
@@ -120,7 +120,7 @@ CreateController(sv, "DamRelease", Y = c("54095", "54057"), U = c("Dam"), FUN =
 
 # GRiwrmRunOptions object simulation of the hydrological year 2002-2003
 IndPeriod_Run <- seq.int(
-  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15")),
+  which(IM_severn[[1]]$DatesR == as.POSIXct("2002-10-15", tz = "UTC")),
   length.out = 365
 )
 IndPeriod_WarmUp <- seq.int(IndPeriod_Run[1] - 366, IndPeriod_Run[1] - 1)
diff --git a/man/plot.GRiwrmOutputsModel.Rd b/man/plot.GRiwrmOutputsModel.Rd
index 60162ec0d3652076fb6744a7b3df8e67857e934e..a2f12a09d49e28228f10befe657c404c61ca4fa0 100644
--- a/man/plot.GRiwrmOutputsModel.Rd
+++ b/man/plot.GRiwrmOutputsModel.Rd
@@ -216,8 +216,8 @@ Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
 
 # Plot the diverted flow for the year 2003
 Ind_Plot <- which(
-  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01") &
-  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31")
+  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
diff --git a/man/plot.Qm3s.Rd b/man/plot.Qm3s.Rd
index acd5b4fa4e04bf30543c0783f48b7bacc4ffbad5..2bb899b633fa59113db059f8936f716c3e485eb2 100644
--- a/man/plot.Qm3s.Rd
+++ b/man/plot.Qm3s.Rd
@@ -247,8 +247,8 @@ Qdiv_m3s <- OM_div$`54001`$Qdiv_m3 / 86400
 
 # Plot the diverted flow for the year 2003
 Ind_Plot <- which(
-  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01") &
-  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31")
+  OM_div[[1]]$DatesR >= as.POSIXct("2003-01-01", tz = "UTC") &
+  OM_div[[1]]$DatesR <= as.POSIXct("2003-12-31", tz = "UTC")
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])