Skip to content

GitLab

Commit a891dd5f authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Accélération d'airGR (modifications des codes fortran des modèles, sauf GR4H) #3832 


Signed-off-by: Delaigue Olivier's avatarOlivier Delaigue <olivier.delaigue@irstea.fr>
parent 21ef1a3f
Hide whitespace changes
Inline Side-by-side
src/frun_GR2M.f
View file @ a891dd5f
......@@ -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
......
src/frun_GR4J.f
View file @ a891dd5f
......@@ -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
......
src/frun_GR5J.f
View file @ a891dd5f
......@@ -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
......
src/frun_GR6J.f
View file @ a891dd5f
......@@ -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
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment