

C**********************************************************************
      SUBROUTINE HU1(XV,C,D)
C Computation of ordinates of GR unit hydrograph HU1 using successives differences on the S curve SS1
C Inputs:
C    C: time constant
C    D: exponent
C Outputs:
C    XV(3*NPX+1) to XV(3*NPX+NH): NH ordinates of discrete hydrograph
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION XV(3*NPX+5*NH)
      DOUBLEPRECISION C,D,SS1
      INTEGER I

      DO I=1,NH
      XV(3*NPX+I)=SS1(I,C,D)-SS1(I-1,C,D)
      ENDDO
      ENDSUBROUTINE


C**********************************************************************
      SUBROUTINE HU2(XV,C,D)
C Computation of ordinates of GR unit hydrograph HU2 using successives differences on the S curve SS2
C Inputs:
C    C: time constant
C    D: exponent
C Outputs:
C    XV(3*NPX+NH+1) to XV(3*NPX+3*NH): 2*NH ordinates of discrete hydrograph
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION XV(3*NPX+5*NH)
      DOUBLEPRECISION C,D,SS2
      INTEGER I

      DO I =1,2*NH
      XV(3*NPX+NH+I)=SS2(I,C,D)-SS2(I-1,C,D)
      ENDDO
      ENDSUBROUTINE



C**********************************************************************
      SUBROUTINE HU4(XV,ALPHA,BETA)
C Computation of ordinates of MOHYSE unit hydrograph
C Inputs:
C    Alpha: parameter
C    Beta:  parameter
C Outputs:
C    XV(3*NPX+NH+1) to XV(3*NPX+3*NH): 2*NH ordinates of discrete hydrograph
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION XV(3*NPX+5*NH),U(3*NH)
      DOUBLEPRECISION ALPHA,BETA,SU
      INTEGER K

      SU=0.
c      IF(ALPHA.LT.1.)THEN  
c      WRITE(*,*)' Pb ALPHA'
c      STOP                 
c      ENDIF                
      IF(ALPHA.EQ.1.)THEN
      U(1)=1.
      SU=1.
      DO 1 K=2,3*NH
      U(K)=0.
    1 CONTINUE
      ELSE
      DO 11 K=1,3*NH
      U(K)=FLOAT(K)*(ALPHA-1.)*EXP(-FLOAT(K)/BETA)
      SU=SU+U(K)
   11 CONTINUE
      ENDIF

c      IF(SU.LT.0.0000000001)THEN       
c      WRITE(*,*)' Pb HU4',ALPHA, BETA  
c      STOP                             
c      ENDIF                            
      DO 2 K=1,3*NH
      XV(3*NPX+K)=U(K)/SU
    2 CONTINUE
      ENDSUBROUTINE



C**********************************************************************
      SUBROUTINE HU(XV,C)
C Computation of ordinates of GRP unit hydrograph
C Inputs:
C    C: time constant
C    Alpha: parameter
C    Beta:  parameter
C Outputs:
C    XV(3*NPX+NH+1) to XV(3*NPX+3*NH): 2*NH ordinates of discrete hydrograph
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION XV(3*NPX+5*NH)
      DOUBLEPRECISION C
      DOUBLEPRECISION SH
      INTEGER I
      DO 10 I=1,2*NH
      XV(3*NPX+NH+I)=SH(I,C)-SH(I-1,C)
   10 CONTINUE
      RETURN
      ENDSUBROUTINE



C**********************************************************************
      FUNCTION SH(I,C)
C Values of the S curve (cumulative HU curve) of GRP unit hydrograph HU
C Inputs:
C    C: time constant
C    I: time-step
C Outputs:
C    SH: Values of the S curve for I
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION C
      DOUBLEPRECISION SH,FI
      INTEGER I
      
      FI=I
      IF(FI.LE.0.)THEN
      SH=0.
      RETURN
      ENDIF
      IF(FI.GE.C)THEN
      SH=1.
      RETURN
      ENDIF
      SH=FI**2.5/(FI**2.5+(C-FI)**2.5)
      RETURN
      ENDFUNCTION


C**********************************************************************
      FUNCTION SS1(I,C,D)
C Values of the S curve (cumulative HU curve) of GR unit hydrograph HU1
C Inputs:
C    C: time constant
C    D: exponent
C    I: time-step
C Outputs:
C    SS1: Values of the S curve for I
C**********************************************************************
      Implicit None
      DOUBLEPRECISION C,D,SS1
      INTEGER I,FI

      FI=I
      IF(FI.LE.0.) THEN
      SS1=0.
      RETURN
      ENDIF
      IF(FI.LT.C) THEN
      SS1=(FI/C)**D
      RETURN
      ENDIF
      SS1=1.
      ENDFUNCTION


C**********************************************************************
      FUNCTION SS2(I,C,D)
C Values of the S curve (cumulative HU curve) of GR unit hydrograph HU2
C Inputs:
C    C: time constant
C    D: exponent
C    I: time-step
C Outputs:
C    SS2: Values of the S curve for I
C**********************************************************************
      Implicit None
      DOUBLEPRECISION C,D,SS2
      INTEGER I,FI

      FI=I
      IF(FI.LE.0.) THEN
      SS2=0.
      RETURN
      ENDIF
      IF(FI.LE.C) THEN
      SS2=0.5*(FI/C)**D
      RETURN
      ENDIF
      IF(FI.LT.2.*C) THEN
      SS2=1.-0.5*(2.-FI/C)**D
      RETURN
      ENDIF
      SS2=1.
      ENDFUNCTION



C**********************************************************************
      SUBROUTINE DEL(XV,C)
C Computation of HU ordinates corresponding to a time lag of a given number (possibly non-integer) of time-steps
C (all ordinates are nul except 2 at max)
C Inputs:
C    C: time constant
C Outputs:
C    XV(3*NPX+NH+1) to XV(3*NPX+3*NH): 2*NH ordinates of discrete hydrograph
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION XV(3*NPX+5*NH)
      DOUBLEPRECISION C,F
      INTEGER I,K
      I=INT(C)
      F=C-INT(C)
      DO 1 K=3*NPX+1,3*NPX+3*NH
      XV(K)=0.
    1 CONTINUE
      XV(3*NPX+I)=1.-F
      XV(3*NPX+I+1)=F
      ENDSUBROUTINE



C**********************************************************************
      SUBROUTINE DEL2(XV,C)
C Computation of HU ordinates corresponding to a time lag of a given number (possibly non-integer) of time-steps
C (all ordinates are nul except 2 at max)
C Inputs:
C    C: time constant
C Outputs:
C    XV(3*NPX+NH+1) to XV(3*NPX+3*NH): NH ordinates of discrete hydrograph
C**********************************************************************
      Implicit None
      INTEGER NPX,NH
      PARAMETER (NPX=14,NH=20)
      DOUBLEPRECISION XV(3*NPX+5*NH)
      DOUBLEPRECISION C,F
      INTEGER K,I
      
      IF(C.GT.FLOAT(NH)) C=FLOAT(NH)
      I=INT(C)
      F=C-INT(C)
      DO 1 K=3*NPX+1,3*NPX+NH
      XV(K)=0.
    1 CONTINUE
      XV(3*NPX+I)=1.-F
      XV(3*NPX+I+1)=F
      ENDSUBROUTINE



C**********************************************************************
      FUNCTION tanHyp(Val)
C Computation of hyperbolic tangent
C**********************************************************************
      Implicit None
      DOUBLEPRECISION Val,ValExp,tanHyp

      ValExp=EXP(Val)
      tanHyp=(ValExp - 1./ValExp)/(ValExp + 1./ValExp)
      RETURN
      ENDFUNCTION

