-
Delaigue Olivier authored76600c38
SUBROUTINE frun_GR6J( !inputs & LInputs , ! [integer] length of input and output series & InputsPrecip , ! [double] input series of total precipitation [mm] & InputsPE , ! [double] input series PE [mm] & NParam , ! [integer] number of model parameter & Param , ! [double] parameter set & NStates , ! [integer] number of state variables used for model initialising & StateStart , ! [double] state variables used when the model run starts (reservoir levels [mm] and HU) & NOutputs , ! [integer] number of output series & IndOutputs , ! [integer] indices of output series !outputs & Outputs , ! [double] output series & StateEnd ) ! [double] state variables at the end of the model run (reservoir levels [mm] and HU) !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr6j Implicit None !### input and output variables integer, intent(in) :: LInputs,NParam,NStates,NOutputs doubleprecision, dimension(LInputs) :: InputsPrecip doubleprecision, dimension(LInputs) :: InputsPE doubleprecision, dimension(NParam) :: Param doubleprecision, dimension(NStates) :: StateStart doubleprecision, dimension(NStates) :: StateEnd integer, dimension(NOutputs) :: IndOutputs doubleprecision, dimension(LInputs,NOutputs) :: Outputs !parameters, internal states and variables integer NPX,NH,NMISC parameter (NPX=14,NH=20,NMISC=30) doubleprecision X(5*NH+7),XV(3*NPX+5*NH) doubleprecision MISC(NMISC) doubleprecision D doubleprecision P1,E,Q integer I,K !-------------------------------------------------------------- !Initialisations !-------------------------------------------------------------- !initilisation of model states to zero X=0. XV=0. !initilisation of model states using StateStart DO I=1,3*NH X(I)=StateStart(I) ENDDO !parameter values !Param(1) : production store capacity (X1 - PROD) [mm] !Param(2) : intercatchment exchange constant (X2 - CES1) [mm/d] !Param(3) : routing store capacity (X3 - ROUT) [mm] !Param(4) : time constant of unit hydrograph (X4 - TB) [d] !Param(5) : intercatchment exchange constant (X5 - CES2) [-] !Param(6) : time constant of exponential store (X6 - EXP) [d] !computation of HU ordinates D=2.5 CALL HU1(XV,Param(4),D) CALL HU2(XV,Param(4),D) !initialisation of model outputs Q = -999.999 MISC = -999.9997172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
c StateEnd = -999.999 !initialisation made in R
c Outputs = -999.999 !initialisation made in R
!--------------------------------------------------------------
!Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P1=InputsPrecip(k)
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
!model run on one time-step
CALL MOD_GR6J(X,XV,Param,P1,E,Q,MISC)
!storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
!model states at the end of the run
DO K=1,3*NH
StateEnd(K)=X(K)
ENDDO
RETURN
ENDSUBROUTINE
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR6J(X,XV,Param,P1,E,Q,MISC)
C Run on a single time-step with the GR6J model
C Inputs:
C X Vector of model states at the beginning of the time-step [mm]
C XV Vector of model states at the beginning of the time-step [mm]
C Param Vector of model parameters [mixed units]
C P1 Value of rainfall during the time-step [mm]
C E Value of potential evapotranspiration during the time-step [mm]
C Outputs:
C X Vector of model states at the end of the time-step [mm]
C XV Vector of model states at the end of the time-step [mm]
C Q Value of simulated flow at the catchment outlet for the time-step [mm]
C MISC Vector of model outputs for the time-step [mm]
C**********************************************************************
Implicit None
INTEGER NPX,NH,NMISC,NParam
PARAMETER (NPX=14,NH=20,NMISC=30)
PARAMETER (NParam=6)
DOUBLEPRECISION X(5*NH+7),XV(3*NPX+5*NH)
DOUBLEPRECISION Param(NParam)
DOUBLEPRECISION MISC(NMISC)
DOUBLEPRECISION P1,E,Q
DOUBLEPRECISION A,B,C,EN,ER,PN,PR,PS,WS,tanHyp,AR
DOUBLEPRECISION PERC,PRHU1,PRHU2,EXCH,QR,QD,QR1
DOUBLEPRECISION AE,AEXCH1,AEXCH2
INTEGER K
DATA B/0.9/
DATA C/0.4/
A=Param(1)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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.
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
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))
X(2)=X(2)-PERC
PR=PR+PERC
PRHU1=PR*B
PRHU2=PR*(1.-B)
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
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
ENDDO
X(7+3*NH)=XV(3*NPX+3*NH)*PRHU2
C Potential intercatchment semi-exchange
EXCH=Param(2)*(X(1)/Param(3)-Param(5))
C Routing store
AEXCH1=EXCH
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.))
X(1)=X(1)-QR
C Update of exponential store
X(6)=X(6)+C*X(8)+EXCH
AR=X(6)/Param(6)
IF(AR.GT.33.)AR=33.
IF(AR.LT.-33.)AR=-33.
IF(AR.GT.7.)THEN
QR1=X(6)+Param(6)/EXP(AR)
GOTO 3
ENDIF
IF(AR.LT.-7.)THEN
QR1=Param(6)*EXP(AR)
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250
GOTO 3
ENDIF
QR1=Param(6)*LOG(EXP(AR)+1.)
3 CONTINUE
X(6)=X(6)-QR1
C Runoff from direct branch QD
AEXCH2=EXCH
IF((X(8+NH)+EXCH).LT.0) AEXCH2=-X(8+NH)
QD=MAX(0.,X(8+NH)+EXCH)
C Total runoff
Q=QR+QD+QR1
IF(Q.LT.0.) Q=0.
C Variables storage
MISC( 1)=E ! PE ! potential evapotranspiration [mm/d]
MISC( 2)=P1 ! Precip ! total precipitation [mm/d]
MISC( 3)=X(2) ! Prod ! production store level (X(2)) [mm]
MISC( 4)=AE ! AE ! actual evapotranspiration [mm/d]
MISC( 5)=PERC ! Perc ! percolation (PERC) [mm]
MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm]
MISC( 7)=X(8) ! Q9 ! outflow from HU1 (Q9) [mm/d]
MISC( 8)=X(8+NH) ! Q1 ! outflow from HU2 (Q1) [mm/d]
MISC( 9)=X(1) ! Rout ! routing store level (X(1)) [mm]
MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/d]
MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/d]
MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/d]
MISC(13)=QR1 ! QR1 ! outflow from exponential store (QR1) [mm/d]
MISC(14)=X(6) ! Exp ! exponential store level (X(6)) (negative) [mm]
MISC(15)=QD ! QD ! outflow from HU2 branch after exchange (QD) [mm/d]
MISC(16)=Q ! Qsim ! outflow at catchment outlet [mm/d]
ENDSUBROUTINE