diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index d61fc315eac94c0ba3d6d1dbfcc39de955050c07..ced3cc714e71f7853fd928799b0ea25a010c7a04 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -132,11 +132,14 @@ CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod
       if(identical(FUN_MOD,RunModel_GR4H)){
         Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
       if(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR4J)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
+        Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                         "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim"); }
       if(identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
+        Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                         "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim"); }
       if(identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
-        Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QR1","Exp","QD","Qsim"); }
+        Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                         "AExch1", "AExch2", "AExch", "QR", "QR1", "Exp", "QD", "Qsim"); }
       if(identical(FUN_MOD,RunModel_GR2M)){
         Outputs_all <- c(Outputs_all,"PotEvap","Precip","AE","Perc","P3","Exch","Prod","Rout","Qsim"); }
       if(identical(FUN_MOD,RunModel_GR1A)){
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index 380225ac050eaae5a91b15991cbeaeff2ca79fe0..1413f96acdd72b707118bdec654fffdc9e5586d4 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -2,7 +2,8 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 6;
     FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
-    FortranOutputsMod       <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
+    FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                        "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index 4a7ea9123e9cf7eab6cd5aa507fcf94d266000b5..a703657b3f698970ecc6c4ef2731cda6e9fec782 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -2,7 +2,8 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 7;
     FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
-    FortranOutputsMod       <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
+    FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                        "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index c4926c650e34e6297923153c140938a3af5b420a..4cf4f11724b35d662b2e5ffddb30efc3f388c4b1 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -2,7 +2,8 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 8;
     FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
-    FortranOutputsMod       <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QR1","Exp","QD","Qsim");
+    FortranOutputsMod       <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1",
+			"Rout", "Exch", "AExch1", "AExch2", "AExch", "QR", "QR1", "Exp", "QD", "Qsim");
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index 24110687ee1651277680403ccb725de3f911af0f..4577fa851b15cf87399307370c224423d8d45c0e 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -1,7 +1,8 @@
 RunModel_GR4J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 4;
-    FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
+    FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                        "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index 3dc9c4630c15d6770da193c127da201d4ffb0d54..6cbf01ff5c36593b70620cc0f3f4d936058f7bff 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -1,7 +1,8 @@
 RunModel_GR5J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 5;
-    FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
+    FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
+                        "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index 9aa0e1b3b1a847b39686dc9aaa711e0ffe9bb676..b7582eb8a231dd586a08fedd9400c9f58ed95c69 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -1,7 +1,8 @@
 RunModel_GR6J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 6;
-    FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QR1","Exp","QD","Qsim");
+    FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1",
+			"Rout", "Exch", "AExch1", "AExch2", "AExch", "QR", "QR1", "Exp", "QD", "Qsim");
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/man/RunModel_CemaNeigeGR4J.Rd b/man/RunModel_CemaNeigeGR4J.Rd
index 209668421518ddf4c92d401837586a1ce9e8a63b..38bbb8518a9cee18c43330341a96411deccf89a8 100644
--- a/man/RunModel_CemaNeigeGR4J.Rd
+++ b/man/RunModel_CemaNeigeGR4J.Rd
@@ -27,6 +27,8 @@ CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d]                       \
          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
          \emph{$Prod    }          \tab [numeric] series of production store level [mm]                        \cr
+         \emph{$Pn      }          \tab [numeric] series of net rainfall [mm/d]                         			  \cr
+         \emph{$Ps      }          \tab [numeric] series of the part of Ps filling the production store [mm/d]        \cr
          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
@@ -34,6 +36,8 @@ CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d]                       \
          \emph{$Q1      }          \tab [numeric] series of UH2 outflow (Q1) [mm/d]                                   \cr
          \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                           \cr
          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
+         \emph{$AExch1  }          \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d]    \cr
+         \emph{$AExch2  }          \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d]    \cr
          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
          \emph{$QD      }          \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d]           \cr
diff --git a/man/RunModel_CemaNeigeGR5J.Rd b/man/RunModel_CemaNeigeGR5J.Rd
index 28fc22195f0fe5c369079ec2e01ceb89b147f82d..01c53c7eb7fcab2104541e70c5e4151c506d70bd 100644
--- a/man/RunModel_CemaNeigeGR5J.Rd
+++ b/man/RunModel_CemaNeigeGR5J.Rd
@@ -28,6 +28,8 @@ CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d]                       \
          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
          \emph{$Prod    }          \tab [numeric] series of production store level [mm]                        \cr
+         \emph{$Pn      }          \tab [numeric] series of net rainfall [mm/d]                         			  \cr
+         \emph{$Ps      }          \tab [numeric] series of the part of Ps filling the production store [mm/d]        \cr
          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
@@ -35,6 +37,8 @@ CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d]                       \
          \emph{$Q1      }          \tab [numeric] series of UH2 outflow (Q1) [mm/d]                                   \cr
          \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                           \cr
          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
+         \emph{$AExch1  }          \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d]    \cr
+         \emph{$AExch2  }          \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d]    \cr
          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
          \emph{$QD      }          \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d]           \cr
diff --git a/man/RunModel_CemaNeigeGR6J.Rd b/man/RunModel_CemaNeigeGR6J.Rd
index cd6e80ff1b39d53881681fc28879e5fcf97ae301..77044b7eb25ebab34eaae9c84605aa60b7cf68b9 100644
--- a/man/RunModel_CemaNeigeGR6J.Rd
+++ b/man/RunModel_CemaNeigeGR6J.Rd
@@ -29,6 +29,8 @@ CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d]                       \
          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
          \emph{$Prod    }          \tab [numeric] series of production store level [mm]                        \cr
+         \emph{$Pn      }          \tab [numeric] series of net rainfall [mm/d]                         			  \cr
+         \emph{$Ps      }          \tab [numeric] series of the part of Ps filling the production store [mm/d]        \cr
          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
@@ -36,6 +38,8 @@ CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d]                       \
          \emph{$Q1      }          \tab [numeric] series of UH2 outflow (Q1) [mm/d]                                   \cr
          \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                           \cr
          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
+         \emph{$AExch1  }          \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d]    \cr
+         \emph{$AExch2  }          \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d]    \cr
          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
          \emph{$QR1     }          \tab [numeric] series of exponential store outflow (QR1) [mm/d]                    \cr
diff --git a/man/RunModel_GR4J.Rd b/man/RunModel_GR4J.Rd
index 89b826df61870f356a29b8b42984e251fb9f4870..c3f488057ae52836ae5e030e6e226219aa8d5522 100644
--- a/man/RunModel_GR4J.Rd
+++ b/man/RunModel_GR4J.Rd
@@ -24,14 +24,18 @@ GR4J X4      \tab unit hydrograph time constant [d]
          \emph{$DatesR  }          \tab [POSIXlt] series of dates                                                     \cr
          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
-         \emph{$Prod    }          \tab [numeric] series of production store level [mm]                        \cr
+         \emph{$Prod    }          \tab [numeric] series of production store level [mm]                 	          \cr
+         \emph{$Pn      }          \tab [numeric] series of net rainfall [mm/d]                         			  \cr
+         \emph{$Ps      }          \tab [numeric] series of the part of Ps filling the production store [mm/d]        \cr
          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
          \emph{$Q9      }          \tab [numeric] series of UH1 outflow (Q9) [mm/d]                                   \cr
          \emph{$Q1      }          \tab [numeric] series of UH2 outflow (Q1) [mm/d]                                   \cr
-         \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                           \cr
+         \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                     		      \cr
          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
+         \emph{$AExch1  }          \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d]    \cr
+         \emph{$AExch2  }          \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d]    \cr
          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
          \emph{$QD      }          \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d]           \cr
diff --git a/man/RunModel_GR5J.Rd b/man/RunModel_GR5J.Rd
index 7c73ba814e846b7563baec4b1ee3468030d3b6fe..803eb8696f76ff274d7c2001b2f03662b51129c4 100644
--- a/man/RunModel_GR5J.Rd
+++ b/man/RunModel_GR5J.Rd
@@ -26,6 +26,8 @@ GR5J X5      \tab intercatchment exchange threshold [-]                     \cr
          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
          \emph{$Prod    }          \tab [numeric] series of production store level [mm]                        \cr
+         \emph{$Pn      }          \tab [numeric] series of net rainfall [mm/d]                         			  \cr
+         \emph{$Ps      }          \tab [numeric] series of the part of Ps filling the production store [mm/d]        \cr
          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
@@ -33,6 +35,8 @@ GR5J X5      \tab intercatchment exchange threshold [-]                     \cr
          \emph{$Q1      }          \tab [numeric] series of UH2 outflow (Q1) [mm/d]                                   \cr
          \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                           \cr
          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
+         \emph{$AExch1  }          \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d]    \cr
+         \emph{$AExch2  }          \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d]    \cr
          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
          \emph{$QD      }          \tab [numeric] series of direct flow from UH2 after exchange (QD) [mm/d]           \cr
diff --git a/man/RunModel_GR6J.Rd b/man/RunModel_GR6J.Rd
index 0727bb5746d0b96dfbafd0a1cd6df8f013f07bd9..5586b906422c2283f87aefeb1ac4b16bf5be75ec 100644
--- a/man/RunModel_GR6J.Rd
+++ b/man/RunModel_GR6J.Rd
@@ -27,6 +27,8 @@ GR6J X6      \tab coefficient for emptying exponential store [mm]
          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
          \emph{$Prod    }          \tab [numeric] series of production store level [mm]                        \cr
+         \emph{$Pn      }          \tab [numeric] series of net rainfall [mm/d]                         			  \cr
+         \emph{$Ps      }          \tab [numeric] series of the part of Ps filling the production store [mm/d]        \cr
          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
@@ -34,6 +36,8 @@ GR6J X6      \tab coefficient for emptying exponential store [mm]
          \emph{$Q1      }          \tab [numeric] series of UH2 outflow (Q1) [mm/d]                                   \cr
          \emph{$Rout    }          \tab [numeric] series of routing store level [mm]                           \cr
          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
+         \emph{$AExch1  }          \tab [numeric] series of actual exchange between catchments for branch 1 [mm/d]    \cr
+         \emph{$AExch2  }          \tab [numeric] series of actual exchange between catchments for branch 2 [mm/d]    \cr
          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
          \emph{$QR1     }          \tab [numeric] series of exponential store outflow (QR1) [mm/d]                    \cr
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index e021eac4a39db30a4773e09d9e1254e4ffd3f2fa..14339005e87c34950b814f745677df0a9b768433 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -173,6 +173,7 @@ C Interception and production store
 	  ! fin speed-up  
       AE=ER+P1
       St(1)=St(1)-ER
+	  PS=0.
       PR=0.
       ELSE
       EN=0.
@@ -254,17 +255,21 @@ C Variables storage
       MISC( 1)=E             ! PE     ! observed potential evapotranspiration [mm/day]
       MISC( 2)=P1            ! Precip ! observed total precipitation [mm/day]
       MISC( 3)=St(1)         ! Prod   ! production store level (St(1)) [mm]
-      MISC( 4)=AE            ! AE     ! actual evapotranspiration [mm/day]
-      MISC( 5)=PERC          ! Perc   ! percolation (PERC) [mm/day]
-      MISC( 6)=PR            ! PR     ! PR=PN-PS+PERC [mm/day]
-      MISC( 7)=StUH1(1)      ! Q9     ! outflow from UH1 (Q9) [mm/day]
-      MISC( 8)=StUH2(1)      ! Q1     ! outflow from UH2 (Q1) [mm/day]
-      MISC( 9)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
-      MISC(10)=EXCH          ! Exch   ! potential semi-exchange between catchments (EXCH) [mm/day]
-      MISC(11)=AEXCH1+AEXCH2 ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
-      MISC(12)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
-      MISC(13)=QD            ! QD     ! outflow from UH2 branch after exchange (QD) [mm/day]
-      MISC(14)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
+      MISC( 4)=PN            ! Pn     ! net rainfall [mm/day]
+      MISC( 5)=PS            ! Ps     ! part of Ps filling the production store [mm/day]
+      MISC( 6)=AE            ! AE     ! actual evapotranspiration [mm/day]
+      MISC( 7)=PERC          ! Perc   ! percolation (PERC) [mm/day]
+      MISC( 8)=PR            ! PR     ! PR=PN-PS+PERC [mm/day]
+      MISC( 9)=StUH1(1)      ! Q9     ! outflow from UH1 (Q9) [mm/day]
+      MISC(10)=StUH2(1)      ! Q1     ! outflow from UH2 (Q1) [mm/day]
+      MISC(11)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
+      MISC(12)=EXCH          ! Exch   ! potential semi-exchange between catchments (EXCH) [mm/day]
+      MISC(13)=AEXCH1 		 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
+      MISC(14)=AEXCH2        ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
+      MISC(15)=AEXCH1+AEXCH2 ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+      MISC(16)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
+      MISC(17)=QD            ! QD     ! outflow from UH2 branch after exchange (QD) [mm/day]
+      MISC(18)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
 
 
 
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index 6744a80bb16c24a45aaa823d9f7f83b62dbe6081..49449d137d5ee5bf2ee8f6f7ae8cf7775a083900 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -163,6 +163,7 @@ C Interception and production store
 	  ! fin speed-up  
       AE=ER+P1
       St(1)=St(1)-ER
+	  PS=0.
       PR=0.
       ELSE
       EN=0.
@@ -239,19 +240,21 @@ C Variables storage
       MISC( 1)=E             ! PE     ! observed potential evapotranspiration [mm/day]
       MISC( 2)=P1            ! Precip ! observed total precipitation [mm/day]
       MISC( 3)=St(1)         ! Prod   ! production store level (St(1)) [mm]
-      MISC( 4)=AE            ! AE     ! actual evapotranspiration [mm/day]
-      MISC( 5)=PERC          ! Perc   ! percolation (PERC) [mm/day]
-      MISC( 6)=PR            ! PR     ! PR=PN-PS+PERC [mm/day]
-      MISC( 7)=Q9            ! Q9     ! outflow from first branch (Q9) [mm/day]
-      MISC( 8)=Q1            ! Q1     ! outflow from second branch (Q1) [mm/day]
-      MISC( 9)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
-      MISC(10)=EXCH          ! Exch   ! potential semi-exchange between catchments (EXCH) [mm/day]
-      MISC(11)=AEXCH1+AEXCH2 ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
-      MISC(12)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
-      MISC(13)=QD            ! QD     ! outflow from UH2 branch after exchange (QD) [mm/day]
-      MISC(14)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
-
-
+      MISC( 4)=PN            ! Pn     ! net rainfall [mm/day]
+      MISC( 5)=PS            ! Ps     ! part of Ps filling the production store [mm/day]
+      MISC( 6)=AE            ! AE     ! actual evapotranspiration [mm/day]
+      MISC( 7)=PERC          ! Perc   ! percolation (PERC) [mm/day]
+      MISC( 8)=PR            ! PR     ! PR=PN-PS+PERC [mm/day]
+      MISC( 9)=Q9	         ! Q9     ! outflow from UH1 (Q9) [mm/day]
+      MISC(10)=Q1		     ! Q1     ! outflow from UH2 (Q1) [mm/day]
+      MISC(11)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
+      MISC(12)=EXCH          ! Exch   ! potential semi-exchange between catchments (EXCH) [mm/day]
+      MISC(13)=AEXCH1 		 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
+      MISC(14)=AEXCH2        ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
+      MISC(15)=AEXCH1+AEXCH2 ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+      MISC(16)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
+      MISC(17)=QD            ! QD     ! outflow from UH2 branch after exchange (QD) [mm/day]
+      MISC(18)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
 
 
       ENDSUBROUTINE
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index 76ca17d6aa1cc661b3d47b7438d13443248401e5..40e316d9a50be55b43f48ed7f8f915723ec5965a 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -180,6 +180,7 @@ C Production store
 	  
       AE=ER+P1
       St(1)=St(1)-ER
+	  PS=0.
       PR=0.
       ELSE
       EN=0.
@@ -283,19 +284,23 @@ C Variables storage
       MISC( 1)=E             ! PE     ! observed potential evapotranspiration [mm/day]
       MISC( 2)=P1            ! Precip ! observed total precipitation [mm/day]
       MISC( 3)=St(1)         ! Prod   ! production store level (St(1)) [mm]
-      MISC( 4)=AE            ! AE     ! actual evapotranspiration [mm/day]
-      MISC( 5)=PERC          ! Perc   ! percolation (PERC) [mm/day]
-      MISC( 6)=PR            ! PR     ! PR=PN-PS+PERC [mm/day]
-      MISC( 7)=StUH1(1)      ! Q9     ! outflow from UH1 (Q9) [mm/day]
-      MISC( 8)=StUH2(1)      ! Q1     ! outflow from UH2 (Q1) [mm/day]
-      MISC( 9)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
-      MISC(10)=EXCH          ! Exch   ! potential third-exchange between catchments (EXCH) [mm/day]
-      MISC(11)=AEXCH1+AEXCH2+EXCH ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2+EXCH) [mm/day]
-      MISC(12)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
-      MISC(13)=QR1           ! QR1    ! outflow from exponential store (QR1) [mm/day]
-      MISC(14)=St(3)         ! Exp    ! exponential store level (St(3)) (negative) [mm]
-      MISC(15)=QD            ! QD     ! outflow from UH2 branch after exchange (QD) [mm/day]
-      MISC(16)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
+      MISC( 4)=PN            ! Pn     ! net rainfall [mm/day]
+      MISC( 5)=PS            ! Ps     ! part of Ps filling the production store [mm/day]
+      MISC( 6)=AE            ! AE     ! actual evapotranspiration [mm/day]
+      MISC( 7)=PERC          ! Perc   ! percolation (PERC) [mm/day]
+      MISC( 8)=PR            ! PR     ! PR=PN-PS+PERC [mm/day]
+      MISC( 9)=StUH1(1)      ! Q9     ! outflow from UH1 (Q9) [mm/day]
+      MISC(10)=StUH2(1)      ! Q1     ! outflow from UH2 (Q1) [mm/day]
+      MISC(11)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
+      MISC(12)=EXCH          ! Exch   ! potential third-exchange between catchments (EXCH) [mm/day]
+      MISC(13)=AEXCH1 		 ! AExch1 ! actual exchange between catchments from routing store (AEXCH1) [mm/day]
+      MISC(14)=AEXCH2        ! AExch2 ! actual exchange between catchments from direct branch (after UH2) (AEXCH2) [mm/day]
+      MISC(15)=AEXCH1+AEXCH2+EXCH ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2+EXCH) [mm/day]
+      MISC(16)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
+      MISC(17)=QR1           ! QR1    ! outflow from exponential store (QR1) [mm/day]
+      MISC(18)=St(3)         ! Exp    ! exponential store level (St(3)) (negative) [mm]
+      MISC(19)=QD            ! QD     ! outflow from UH2 branch after exchange (QD) [mm/day]
+      MISC(20)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
 
 
       ENDSUBROUTINE