diff --git a/DESCRIPTION b/DESCRIPTION
index 8b40fa90aa3ac4f463dee0a4998f80071ebe7fe1..023cb6d52a6ecda6617fb2e376635e866943555e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.3.2.64
-Date: 2019-11-21
+Version: 1.3.2.65
+Date: 2019-12-02
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
   person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@irstea.fr"),
diff --git a/NEWS.md b/NEWS.md
index 4f4bc8b15c4b428877ea97c231624c2ed7f7ff1e..d7007eeebadb615e271f494e4e7e0089f66ba109 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 
 
-### 1.3.2.64 Release Notes (2019-11-21)
+### 1.3.2.65 Release Notes (2019-12-02)
 
 
 #### New features
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f
index 2789e0b9cd0304ea267842c1fac915dac4b88af7..e35857c6b043cb19a84e4fa47d7891772824217a 100644
--- a/src/frun_GR1A.f
+++ b/src/frun_GR1A.f
@@ -11,7 +11,7 @@
 ! Further cleaning: G. Thirel
 !------------------------------------------------------------------------------
 ! Creation date: 2003
-! Last modified: 21/11/2019
+! Last modified: 25/11/2019
 !------------------------------------------------------------------------------
 ! REFERENCES
 ! Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit 
@@ -97,7 +97,7 @@
         CALL MOD_GR1A(Param,P0,P1,E1,Q,MISC)
         ! storage of outputs
         DO I=1,NOutputs
-        Outputs(k,I)=MISC(IndOutputs(I))
+          Outputs(k,I)=MISC(IndOutputs(I))
         ENDDO
       ENDDO
 
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index d23dc3b9590ee2b8765cd1fddd09132ca283ffc7..c3229b592e6b9a8320efc8fea16d3c49a7e842a3 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -11,7 +11,7 @@
 ! Further cleaning: G. Thirel
 !------------------------------------------------------------------------------
 ! Creation date: 2003
-! Last modified: 21/11/2019
+! Last modified: 25/11/2019
 !------------------------------------------------------------------------------
 ! REFERENCES
 ! Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit 
@@ -161,7 +161,7 @@
       
 ! Production store
       WS=P/Param(1)  
-      IF(WS.GT.13.)WS=13.
+      IF(WS.GT.13.) WS=13.
 
       ! speed-up
       TWS = tanHyp(WS)
@@ -171,7 +171,7 @@
 
       P1=P+St(1)-S1                  
       WS=E/Param(1)         
-      IF(WS.GT.13.)WS=13.
+      IF(WS.GT.13.) WS=13.
 
       ! speed-up
       TWS = tanHyp(WS)
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f
index 03cf31cce165877c9217631c0ed4083f50796e17..a63ccedc6f368934f3d76f3703d67166ea24e200 100644
--- a/src/frun_GR4H.f
+++ b/src/frun_GR4H.f
@@ -11,7 +11,7 @@
 ! Further cleaning: G. Thirel
 !------------------------------------------------------------------------------
 ! Creation date: 2003
-! Last modified: 21/11/2019
+! Last modified: 25/11/2019
 !------------------------------------------------------------------------------
 ! REFERENCES
 ! Perrin, C., C. Michel and V. Andréassian (2003). Improvement of a 
@@ -200,10 +200,10 @@
 
 ! Interception and production store
       IF(P1.LE.E) THEN
-      EN=E-P1
-      PN=0.
-      WS=EN/A
-      IF(WS.GT.13.)WS=13.
+        EN=E-P1
+        PN=0.
+        WS=EN/A
+        IF(WS.GT.13.) WS=13.
       
       ! speed-up
         TWS = tanHyp(WS)
@@ -212,15 +212,15 @@
       ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
       ! end speed-up  
       
-      AE=ER+P1
-      St(1)=St(1)-ER
-      PR=0.
+        AE=ER+P1
+        St(1)=St(1)-ER
+        PR=0.
       ELSE
-      EN=0.
-      AE=E
-      PN=P1-E
-      WS=PN/A
-      IF(WS.GT.13.)WS=13.
+        EN=0.
+        AE=E
+        PN=P1-E
+        WS=PN/A
+        IF(WS.GT.13.)WS=13.
       
       ! speed-up
         TWS = tanHyp(WS)
@@ -229,19 +229,19 @@
       ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
       ! end speed-up
       
-      PR=PN-PS
-      St(1)=St(1)+PS
+        PR=PN-PS
+        St(1)=St(1)+PS
       ENDIF
 
 ! Percolation from production store
-      IF(St(1).LT.0.)St(1)=0.
+      IF(St(1).LT.0.) St(1)=0.
 
       ! speed-up
-      ! (21/4)**4 = 759.69140625
+      ! (21/4)**4 = 759.69140625 = stored_val
       Sr = St(1)/Param(1)
       Sr = Sr * Sr
       Sr = Sr * Sr
-      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/759.69140625)))
+      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val)))
       ! PERC=X(2)*(1.-(1.+(X(2)/(21./4.*Param(1)))**4.)**(-0.25))
       ! end speed-up
 
@@ -255,13 +255,13 @@
 
 ! Convolution of unit hydrograph UH1
       DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
-      StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
+        StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
       ENDDO
       StUH1(NH)=OrdUH1(NH)*PRHU1
 
 ! Convolution of unit hydrograph UH2
       DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1.)))
-      StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRHU2
+        StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRHU2
       ENDDO
       StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
 
@@ -276,7 +276,7 @@
       AEXCH1=EXCH
       IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
       St(2)=St(2)+StUH1(1)+EXCH
-      IF(St(2).LT.0.)St(2)=0.
+      IF(St(2).LT.0.) St(2)=0.
 
       ! speed-up
       Rr = St(2)/Param(3)
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index c3a88d4da3b906600269d5b52c1a34d28d68b08b..c92837e93071ced1aaa36ca3d442e67acc2287c6 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -11,7 +11,7 @@
 ! Further cleaning: G. Thirel
 !------------------------------------------------------------------------------
 ! Creation date: 2000
-! Last modified: 21/11/2019
+! Last modified: 25/11/2019
 !------------------------------------------------------------------------------
 ! REFERENCES
 ! Perrin, C., C. Michel and V. Andréassian (2003). Improvement of a 
@@ -84,10 +84,10 @@
       St(1) = StateStart(1)
       St(2) = StateStart(2)
       DO I=1,NH
-      StUH1(I)=StateStart(7+I)
+        StUH1(I)=StateStart(7+I)
       ENDDO
       DO I=1,2*NH
-      StUH2(I)=StateStart(7+I+NH)
+        StUH2(I)=StateStart(7+I+NH)
       ENDDO
 
       ! parameter values
@@ -124,17 +124,17 @@
         CALL MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
         ! storage of outputs
         DO I=1,NOutputs
-        Outputs(k,I)=MISC(IndOutputs(I))
+          Outputs(k,I)=MISC(IndOutputs(I))
         ENDDO
       ENDDO
       ! model states at the end of the run
       StateEnd(1) = St(1)
       StateEnd(2) = St(2)
       DO K=1,NH
-      StateEnd(7+K)=StUH1(K)
+        StateEnd(7+K)=StUH1(K)
       ENDDO
       DO K=1,2*NH
-      StateEnd(7+NH+K)=StUH2(K)
+        StateEnd(7+NH+K)=StUH2(K)
       ENDDO
 
       RETURN
@@ -200,44 +200,44 @@
 
 ! Interception and production store
       IF(P1.LE.E) THEN
-      EN=E-P1
-      PN=0.
-      WS=EN/A
-      IF(WS.GT.13.)WS=13.
-      ! speed-up
-      TWS = tanHyp(WS)
-      Sr = St(1)/A
-      ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
-      ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
-      ! end speed-up  
-      AE=ER+P1
-      St(1)=St(1)-ER
-      PS=0.
-      PR=0.
+        EN=E-P1
+        PN=0.
+        WS=EN/A
+        IF(WS.GT.13.) WS=13.
+        ! speed-up
+        TWS = tanHyp(WS)
+        Sr = St(1)/A
+        ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+        ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+        ! end speed-up  
+        AE=ER+P1
+        St(1)=St(1)-ER
+        PS=0.
+        PR=0.
       ELSE
-      EN=0.
-      AE=E
-      PN=P1-E
-      WS=PN/A
-      IF(WS.GT.13.)WS=13.
-      ! speed-up
-      TWS = tanHyp(WS)
-      Sr = St(1)/A
-      PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
-      ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
-      ! end speed-up
-      PR=PN-PS
-      St(1)=St(1)+PS
+        EN=0.
+        AE=E
+        PN=P1-E
+        WS=PN/A
+        IF(WS.GT.13.) WS=13.
+        ! speed-up
+        TWS = tanHyp(WS)
+        Sr = St(1)/A
+        PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
+        ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
+        ! end speed-up
+        PR=PN-PS
+        St(1)=St(1)+PS
       ENDIF
 
 ! Percolation from production store
-      IF(St(1).LT.0.)St(1)=0.
+      IF(St(1).LT.0.) St(1)=0.
       ! speed-up
-      ! (9/4)**4 = 25.62891
+      ! (9/4)**4 = 25.62890625 = stored_val
       Sr = St(1)/Param(1)
       Sr = Sr * Sr
       Sr = Sr * Sr
-      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62891)))
+      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val)))
       ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
       ! end speed-up
       St(1)=St(1)-PERC
@@ -250,13 +250,13 @@
 
 ! Convolution of unit hydrograph UH1
       DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
-      StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
+        StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
       ENDDO
       StUH1(NH)=OrdUH1(NH)*PRHU1
 
 ! Convolution of unit hydrograph UH2
       DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1.)))
-      StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRHU2
+        StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRHU2
       ENDDO
       StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
 
@@ -271,7 +271,7 @@
       AEXCH1=EXCH
       IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
       St(2)=St(2)+StUH1(1)+EXCH
-      IF(St(2).LT.0.)St(2)=0.
+      IF(St(2).LT.0.) St(2)=0.
       ! speed-up
       Rr = St(2)/Param(3)
       Rr = Rr * Rr
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index a5baab6a2c279a63ef7536efab565ee4d73fc2a8..904d31e54c8c5fe32dcc458bd95a1c4d739f707b 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -11,7 +11,7 @@
 ! Further cleaning: G. Thirel
 !------------------------------------------------------------------------------
 ! Creation date: 2006
-! Last modified: 21/11/2019
+! Last modified: 25/11/2019
 !------------------------------------------------------------------------------
 ! REFERENCES
 ! Le Moine, N. (2008). Le bassin versant de surface vu par le souterrain : une 
@@ -88,7 +88,7 @@
       St(2) = StateStart(2)
 
       DO I=1,2*NH
-      StUH2(I)=StateStart(7+NH+I)
+        StUH2(I)=StateStart(7+NH+I)
       ENDDO
 
       ! parameter values
@@ -124,14 +124,14 @@
         CALL MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,MISC)
         ! storage of outputs
         DO I=1,NOutputs
-        Outputs(k,I)=MISC(IndOutputs(I))
+          Outputs(k,I)=MISC(IndOutputs(I))
         ENDDO
       ENDDO
       ! model states at the end of the run
       StateEnd(1) = St(1)
       StateEnd(2) = St(2)
       DO K=1,2*NH
-      StateEnd(7+NH+K)=StUH2(K)
+        StateEnd(7+NH+K)=StUH2(K)
       ENDDO
 
       RETURN
@@ -192,10 +192,10 @@
 
 ! Interception and production store
       IF(P1.LE.E) THEN
-      EN=E-P1
-      PN=0.
-      WS=EN/A
-      IF(WS.GT.13.)WS=13.
+        EN=E-P1
+        PN=0.
+        WS=EN/A
+        IF(WS.GT.13.) WS=13.
         ! speed-up
         TWS = tanHyp(WS)
         Sr = St(1)/A
@@ -211,7 +211,7 @@
         AE=E
         PN=P1-E
         WS=PN/A
-        IF(WS.GT.13.)WS=13.
+        IF(WS.GT.13.) WS=13.
         ! speed-up
         TWS = tanHyp(WS)
         Sr = St(1)/A
@@ -224,14 +224,14 @@
       ENDIF
 
 ! Percolation from production store
-      IF(St(1).LT.0.)St(1)=0.
+      IF(St(1).LT.0.) St(1)=0.
 
       ! speed-up
-      ! (9/4)**4 = 25.62890625
+      ! (9/4)**4 = 25.62890625 = stored_val
       Sr = St(1)/Param(1)
       Sr = Sr * Sr
       Sr = Sr * Sr
-      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val)))
       ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
       ! end speed-up
 
@@ -256,13 +256,13 @@
       AEXCH1=EXCH
       IF((St(2)+Q9+EXCH).LT.0.) AEXCH1=-St(2)-Q9
       St(2)=St(2)+Q9+EXCH
-      IF(St(2).LT.0.)St(2)=0.
+      IF(St(2).LT.0.) St(2)=0.
 
       ! speed-up
       Rr = St(2)/Param(3)
       Rr = Rr * Rr
       Rr = Rr * Rr
-      QR=St(2)*(1.-1./SQRT(SQRT(1.+Rr)))
+      QR = St(2)*(1.-1./SQRT(SQRT(1.+Rr)))
       ! QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
       ! end speed-up
 
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index d2c7d2cdf8b8aa94d2b4fcb6e0e9071212719dd6..5e3bc84a3f761da6b23cd71b56d694db56227083 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -11,7 +11,7 @@
 ! Further cleaning: G. Thirel
 !------------------------------------------------------------------------------
 ! Creation date: 2010
-! Last modified: 21/11/2019
+! Last modified: 25/11/2019
 !------------------------------------------------------------------------------
 ! REFERENCES
 ! Pushpalatha, R., C. Perrin, N. Le Moine, T. Mathevet and V. Andréassian! 
@@ -209,7 +209,7 @@
         EN=E-P1
         PN=0.
         WS=EN/A
-        IF(WS.GT.13)WS=13.
+        IF(WS.GT.13) WS=13.
         
         ! speed-up
         TWS = tanHyp(WS)
@@ -227,7 +227,7 @@
         AE=E
         PN=P1-E
         WS=PN/A
-        IF(WS.GT.13)WS=13.
+        IF(WS.GT.13) WS=13.
         
         ! speed-up
         TWS = tanHyp(WS)
@@ -241,13 +241,13 @@
       ENDIF
 
 ! Percolation from production store
-      IF(St(1).LT.0.)St(1)=0.
+      IF(St(1).LT.0.) St(1)=0.
       ! speed-up
-      ! (9/4)**4 = 25.62890625
+      ! (9/4)**4 = 25.62890625 = stored_val
       Sr = St(1)/Param(1)
       Sr = Sr * Sr
       Sr = Sr * Sr
-      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+      PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val)))
       ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
       ! end speed-up  
 
@@ -261,13 +261,13 @@
 
 ! Convolution of unit hydrograph UH1
       DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
-      StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRUH1
+        StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRUH1
       ENDDO
       StUH1(NH)=OrdUH1(NH)*PRUH1
 
 ! Convolution of unit hydrograph UH2
       DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1.)))
-      StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRUH2
+        StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRUH2
       ENDDO
       StUH2(2*NH)=OrdUH2(2*NH)*PRUH2
 
@@ -278,7 +278,7 @@
       AEXCH1=EXCH
       IF((St(2)+(1-C)*StUH1(1)+EXCH).LT.0) AEXCH1=-St(2)-(1-C)*StUH1(1)
       St(2)=St(2)+(1-C)*StUH1(1)+EXCH
-      IF(St(2).LT.0.)St(2)=0.
+      IF(St(2).LT.0.) St(2)=0.
       
       ! speed-up
       Rr = St(2)/Param(3)
@@ -293,16 +293,16 @@
 ! Update of exponential store
       St(3)=St(3)+C*StUH1(1)+EXCH
       AR=St(3)/Param(6)
-      IF(AR.GT.33.)AR=33.
-      IF(AR.LT.-33.)AR=-33.
+      IF(AR.GT.33.)  AR=33.
+      IF(AR.LT.-33.) AR=-33.
 
       IF(AR.GT.7.)THEN
-      QRExp=St(3)+Param(6)/EXP(AR)
+        QRExp=St(3)+Param(6)/EXP(AR)
       GOTO 3
       ENDIF
 
       IF(AR.LT.-7.)THEN
-      QRExp=Param(6)*EXP(AR)
+        QRExp=Param(6)*EXP(AR)
       GOTO 3
       ENDIF