diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index a7ae201ceb708455f88cb41e8d84dea19ba59f9d..b70245c0bf8c3ea411f21423895a121b723c83b8 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -117,18 +117,36 @@ C**********************************************************************
       DOUBLEPRECISION WS,tanHyp,S1,S2
       DOUBLEPRECISION P1,P2,P3,R1,R2
 
-
+      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
+	  
 C Production store
       WS=P/Param(1)  
-      IF(WS.GT.13)WS=13    
-      S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))                 
-      P1=P+X(1)-S1                  
+      IF(WS.GT.13)WS=13   
+	  
+ 	  ! speed-up
+	  TWS = tanHyp(WS)
+	  S1=(X(1)+Param(1)*TWS)/(1.+X(1)/Param(1)*TWS)                 
+	  ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))                 
+  	  ! fin speed-up
+
+	  P1=P+X(1)-S1                  
       WS=E/Param(1)         
-      IF(WS.GT.13)WS=13    
-      S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS))  
+      IF(WS.GT.13)WS=13  
+	  
+ 	  ! speed-up
+	  TWS = tanHyp(WS)
+      S2=S1*(1.-TWS)/(1.+(1.-S1/Param(1))*TWS)  
+      ! S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS))  
+ 	  ! fin speed-up
                 
 C Percolation
-      X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)         
+ 	  ! speed-up
+	  Sr = S2/Param(1)
+	  Sr = Sr * Sr * Sr + 1.
+      X(1)=S2/Sr**(1./3.)
+      ! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)         
+ 	  ! fin speed-up
+
       P2=S2-X(1)  
       P3=P1+P2
 
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index 953e7647eb4dcf2fed539bf9543f96b4fb30ce50..8db23f835e7052cf5ed61671274d9a189686d2b3 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -133,35 +133,60 @@ C**********************************************************************
       INTEGER K
 
       DATA B/0.9/
+	  
+      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
 
       A=Param(1)
 
 
 C Production store
       IF(P1.LE.E) THEN
-      EN=E-P1
-      PN=0.
-      WS=EN/A
-      IF(WS.GT.13)WS=13.
-      ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
-      AE=ER+P1
-      IF(X(2).LT.ER) AE=X(2)+P1
-      X(2)=X(2)-ER
-      PR=0.
+		EN=E-P1
+		PN=0.
+		WS=EN/A
+		IF(WS.GT.13)WS=13.
+	  
+	  ! speed-up
+		TWS = tanHyp(WS)
+		Sr = X(2)/A
+		ER=X(2)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+      ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+	  ! fin speed-up  
+	  
+		AE=ER+P1
+		IF(X(2).LT.ER) AE=X(2)+P1
+		X(2)=X(2)-ER
+		PR=0.
       ELSE
-      EN=0.
-      AE=E
-      PN=P1-E
-      WS=PN/A
-      IF(WS.GT.13)WS=13.
-      PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
-      PR=PN-PS
-      X(2)=X(2)+PS
+		EN=0.
+		AE=E
+		PN=P1-E
+		WS=PN/A
+		IF(WS.GT.13)WS=13.
+	  
+	  ! speed-up
+		TWS = tanHyp(WS)
+		Sr = X(2)/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))
+	  ! fin speed-up
+	  
+		PR=PN-PS
+		X(2)=X(2)+PS
       ENDIF
 
 C Percolation from production store
       IF(X(2).LT.0.)X(2)=0.
-      PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+
+	  ! speed-up
+	  ! (9/4)**4 = 25.62891
+ 	  Sr = X(2)/Param(1)
+	  Sr = Sr * Sr
+	  Sr = Sr * Sr
+      PERC=X(2)*(1.-1./SQRT(SQRT(1.+Sr/25.62891)))
+	  ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+	  ! fin speed-up
+
       X(2)=X(2)-PERC
 
       PR=PR+PERC
@@ -171,25 +196,37 @@ C Percolation from production store
 
 C Unit hydrograph HU1
       DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1)))
-      X(7+K)=X(8+K)+XV(3*NPX+K)*PRHU1
+		X(7+K)=X(8+K)+XV(3*NPX+K)*PRHU1
       ENDDO
       X(7+NH)=XV(3*NPX+NH)*PRHU1
 
 C Unit hydrograph HU2
       DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1)))
-      X(7+NH+K)=X(8+NH+K)+XV(3*NPX+NH+K)*PRHU2
+		X(7+NH+K)=X(8+NH+K)+XV(3*NPX+NH+K)*PRHU2
       ENDDO
       X(7+3*NH)=XV(3*NPX+3*NH)*PRHU2
 
 C Potential intercatchment semi-exchange
-      EXCH=Param(2)*(X(1)/Param(3))**3.5
+	  ! speed-up
+	  Rr = X(1)/Param(3)
+      EXCH=Param(2)*Rr*Rr*Rr*SQRT(Rr)
+      ! EXCH=Param(2)*(X(1)/Param(3))**3.5
+	  ! fin speed-up
 
 C Routing store
       AEXCH1=EXCH
       IF((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
       X(1)=X(1)+X(8)+EXCH
       IF(X(1).LT.0.)X(1)=0.
-      QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
+
+	  ! speed-up
+	  Rr = X(1)/Param(3)
+	  Rr = Rr * Rr
+	  Rr = Rr * Rr
+      QR=X(1)*(1.-1./SQRT(SQRT(1.+Rr)))
+      ! QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
+	  ! fin speed-up
+
       X(1)=X(1)-QR
 
 C Runoff from direct branch QD
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index 4b7b73515685d0ebeaf68a0b43cfc70d3978b320..b8fb3d02ea6dd9aaf1272f3b08aba5ba0c8ee861 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -134,6 +134,8 @@ C**********************************************************************
       INTEGER K
 
       DATA B/0.9/
+	  
+      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
 
       A=Param(1)
 
@@ -144,7 +146,14 @@ C Production store
       PN=0.
       WS=EN/A
       IF(WS.GT.13)WS=13.
-      ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+	  
+	  ! speed-up
+	  TWS = tanHyp(WS)
+	  Sr = X(2)/A
+      ER=X(2)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+      ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+	  ! fin speed-up  
+	  
       AE=ER+P1
       IF(X(2).LT.ER) AE=X(2)+P1
       X(2)=X(2)-ER
@@ -155,14 +164,30 @@ C Production store
       PN=P1-E
       WS=PN/A
       IF(WS.GT.13)WS=13.
-      PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
+	  
+	  ! speed-up
+	  TWS = tanHyp(WS)
+	  Sr = X(2)/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))
+	  ! fin speed-up
+	  
       PR=PN-PS
       X(2)=X(2)+PS
       ENDIF
 
 C Percolation from production store
       IF(X(2).LT.0.)X(2)=0.
-      PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+
+	  ! speed-up
+	  ! (9/4)**4 = 25.62890625
+ 	  Sr = X(2)/Param(1)
+	  Sr = Sr * Sr
+	  Sr = Sr * Sr
+      PERC=X(2)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+	  ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+	  ! fin speed-up
+
       X(2)=X(2)-PERC
 
       PR=PR+PERC
@@ -190,7 +215,15 @@ C Routing store
       IF((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
       X(1)=X(1)+X(8)+EXCH
       IF(X(1).LT.0.)X(1)=0.
-      QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
+
+	  ! speed-up
+	  Rr = X(1)/Param(3)
+	  Rr = Rr * Rr
+	  Rr = Rr * Rr
+      QR=X(1)*(1.-1./SQRT(SQRT(1.+Rr)))
+      ! QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
+	  ! fin speed-up
+
       X(1)=X(1)-QR
 
 C Runoff from direct branch QD
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index 9d7f9373ffb61bb0850481dd3c92ad132aa999de..06ac8253a28f1a7aeaa1d4659a7ea42c5583cfce 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -136,6 +136,8 @@ C**********************************************************************
 
       DATA B/0.9/
       DATA C/0.4/
+	  
+      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
 
       A=Param(1)
 
@@ -146,7 +148,14 @@ C Production store
       PN=0.
       WS=EN/A
       IF(WS.GT.13)WS=13.
-      ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+	  
+	  ! speed-up
+	  TWS = tanHyp(WS)
+	  Sr = X(2)/A
+      ER=X(2)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+      ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+	  ! fin speed-up  
+	  
       AE=ER+P1
       IF(X(2).LT.ER) AE=X(2)+P1
       X(2)=X(2)-ER
@@ -157,14 +166,29 @@ C Production store
       PN=P1-E
       WS=PN/A
       IF(WS.GT.13)WS=13.
-      PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
+	  
+	  ! speed-up
+	  TWS = tanHyp(WS)
+	  Sr = X(2)/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))
+	  ! fin speed-up  
+	  
       PR=PN-PS
       X(2)=X(2)+PS
       ENDIF
 
 C Percolation from production store
       IF(X(2).LT.0.)X(2)=0.
-      PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+  	  ! speed-up
+	  ! (9/4)**4 = 25.62890625
+ 	  Sr = X(2)/Param(1)
+	  Sr = Sr * Sr
+	  Sr = Sr * Sr
+      PERC=X(2)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+      ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+	  ! fin speed-up  
+
       X(2)=X(2)-PERC
 
       PR=PR+PERC
@@ -192,7 +216,15 @@ C Routing store
       IF((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
       X(1)=X(1)+(1-C)*X(8)+EXCH
       IF(X(1).LT.0.)X(1)=0.
-      QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
+	  
+	  ! speed-up
+	  Rr = X(1)/Param(3)
+	  Rr = Rr * Rr
+	  Rr = Rr * Rr
+      QR=X(1)*(1.-1./SQRT(SQRT(1.+Rr)))
+      ! QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
+	  ! fin speed-up
+	  
       X(1)=X(1)-QR
 
 C Update of exponential store