C
C**************************  MC ****************************************
C  MUSKINGUM CUNGE METHOD: FROM THE INFLOW HYDROGRAPHS AT T+1 AND T,
C  Q_IN_KKP1 AND Q_IN_KK RESPECTIVELY, FROM THE OUTFLOW HYDROGRAPH AT T,
C  Q_OUT_KK, AND FROM THE LATERAL CONTRIBUTION, Q_OVERLAND, IT DETERMINES 
C  THE OUTFLOW HYDROGRAPH FOR THE DOWNSTREAM CELL AT T+1, Q_OUT_KKP1.
C***********************************************************************
C

      SUBROUTINE MC(LOCAL_SLOPE,EPL,KS,W,B1,Y1,DELTAT,CU,PE,X,C1,
     1     AK,Q_IN_KK,Q_IN_KKP1,Q_OUT_KK,Q_OVERLAND,Q_OUT_KKP1)    
              
      IMPLICIT NONE
      
      REAL*8 LOCAL_SLOPE,EPL,KS,W
      REAL*8 B1,Y1,DELTAT
      REAL*8 CU
      REAL*8 Q_IN_KK,Q_IN_KKP1,Q_OUT_KK,Q_OUT_KKP1,Q_OVERLAND
      REAL*8 QC,QCMIN,BETA,G
      REAL*8 CK,DH,PE,X,AK
      REAL*8 C1,C2,C3,C4
      
      
C     THREE-POINT AVERAGE DISCHARGE
      QC=1.0D0/3.0D0*(Q_IN_KK+Q_IN_KKP1+Q_OUT_KK)
      QCMIN=1.0D-05
      IF(QC.LE.QCMIN)  QC=QCMIN
      
      
C     SLOPE
      BETA=ATAN(LOCAL_SLOPE)

C     HYDRAULIC GEOMETRY TERM
      G=(1.0D0-Y1+2.0D0/3.0D0*B1)

C     KINEMATIC WAVE CELERITY
      CK=5.0D0/(3.0D0*G)*KS**(3.0D0/5.0D0)*W**(-2.0D0/5.0D0)*
     1     SIN(BETA)**(3.0D0/1.0D1)*QC**(1.0D0-3.0D0*G/5.0D0)       

      AK=CK/EPL

C     COURANT NUMBER
      CU=CK*DELTAT/EPL

C     HYDRAULIC DIFFUSIVITY
      DH=(QC**(1.0D0-B1))/(2*G*W*TAN(BETA))      
      
      IF (DH.LT.(1.0D0-CU)) DH=1.0D0-CU
      
C     PECLET NUMBER
      PE=CK*EPL/DH

      X=0.50D0-DH/(CK*EPL)

      C1=(CU-2.0D0*X)/(2.0D0*(1.0D0-X)+CU)

      C2=(CU+2.0D0*X)/(2.0D0*(1.0D0-X)+CU)

      C3=(2.0D0*(1.0D0-X)-CU)/(2.0D0*(1.0D0-X)+CU)

      C4=(2.0D0*CK*DELTAT)/(2.0D0*(1.0D0-X)+CU)  

      Q_OUT_KKP1=C1*Q_IN_KKP1+C2*Q_IN_KK+C3*Q_OUT_KK+C4*Q_OVERLAND

      RETURN
      END
