diff --git a/NEWS.md b/NEWS.md
index 722bd0e1cd01e2bb5bde3f63404450c4f2e584c4..bcdf86dd449a5db10564103f38653464065cd1e4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,6 +3,9 @@
 - New features:
     + PlotPDFsOLE includes parameters to modify legend style
 
+- Fixes:
+    + PlotForecastPDF correctly displays terciles labels 
+
 ### CSTools 3.1.0
 **Submission date to CRAN: 02-07-2020**
 
diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R
index 14e6cd3aadca4cda3f6d040e526b8bcc5a0b206a..9cd4dfebcad4c7b037daacd4d681ed534289f55e 100644
--- a/R/PlotForecastPDF.R
+++ b/R/PlotForecastPDF.R
@@ -49,19 +49,19 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
     #------------------------
     color.set <- match.arg(color.set)
     if (color.set == "s2s4e") {
-        colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1"))
+        colorFill <- c("#FF764D", "#b5b5b5", "#33BFD1")
         colorHatch <- c("deepskyblue3", "indianred3")
         colorMember <- c("#ffff7f")
         colorObs <- "purple"
         colorLab <- c("blue", "red")
     } else if (color.set == "hydro") {
-        colorFill <- rev(c("#41CBC9", "#b5b5b5", "#FFAB38"))
+        colorFill <- c("#41CBC9", "#b5b5b5", "#FFAB38")
         colorHatch <- c("darkorange1", "deepskyblue3")
         colorMember <- c("#ffff7f")
         colorObs <- "purple"
         colorLab <- c("darkorange3", "blue")
     } else if (color.set == "ggplot") {
-        colorFill <- rev(ggColorHue(3))
+        colorFill <- ggColorHue(3)
         colorHatch <- c("deepskyblue3", "indianred1")
         colorMember <- c("#ffff7f")
         colorObs <- "purple"
@@ -306,6 +306,12 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
     tmp.dt <- data.table(tmp.df)
     pct <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), 
         by = .(init, tercile)]
+    # include potentially missing groups
+    pct <- merge(pct, CJ(init = factor(levels(pct$init), levels = levels(pct$init)), 
+                         tercile = factor(c("Below normal", "Normal", "Above normal"), 
+                         levels = c("Below normal", "Normal", "Above normal"))),
+                 by = c("init", "tercile"), all.y = T)
+    pct[is.na(pct),"pct"] <- 0
     tot <- pct[, .(tot = sum(pct)), by = init]
     pct <- merge(pct, tot, by = "init")
     pct$pct <- round(100 * pct$pct/pct$tot, 0)
@@ -317,15 +323,17 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
     if (!is.null(extreme.limits)) {
         pct2 <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), 
             by = .(init, extremes)]
+        # include potentially missing groups
+        pct2 <- merge(pct2, CJ(init = factor(levels(pct2$init), levels = levels(pct2$init)), 
+                               extremes = factor(c("Below P10", "Normal", "Above P90"),
+                               levels = c("Below P10", "Normal", "Above P90"))),
+                      by = c("init", "extremes"), all.y=T)
+        pct2[is.na(pct),"pct"] <- 0
         tot2 <- pct2[, .(tot = sum(pct)), by = init]
         pct2 <- merge(pct2, tot2, by = "init")
         pct2$pct <- round(100 * pct2$pct/pct2$tot, 0)
         pct2$lab.pos <-  as.vector(apply(extreme.limits, 1, function(x) {c(x[1], NA, x[2])}))
-        pct2 <- merge(pct2, max.df, by = c("init", "extremes"))
-        # include potentially missing groups
-        pct2 <- pct2[CJ(factor(levels(pct2$init), levels = levels(pct2$init)), 
-            factor(c("Below P10", "Normal", "Above P90"), 
-            levels = c("Below P10", "Normal", "Above P90"))), ]
+        pct2 <- merge(pct2, max.df, by = c("init", "extremes"), all.x = T)
     }
     #------------------------
     # Add probability labels for terciles