

      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.999
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)


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)
      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


