      SUBROUTINE CHANNEL_INC1(DRAIN_INCISE,OUT_I,OUT_J,DEST_I,DEST_J,
     &RUNOFF,LAKE_HT,LAKE_AR,PP,QQ,LAKE_BASE,SURF,DEM,PATH_TOTAL,MCI,
     &CHN_DEC,SUMINC,RUNOFF_X,RUNOFF_Y,RUNOFF_LAK,XIND,YIND,LAKE_RED,K1
     &,K2,K3,K4,CMARK,LAKE2RUN,ROWS,COLS,SX)

      ! THIS SUBROUTINE CALCULATES THE DEPTH OF THE CHANNEL INCISED
      ! DURING LAKE OVERFLOW

      ! DESCRIPTION OF VARIABLES:

      ! DRAIN_INCISE: STORES THE LOCATION OF THE NUCLEI OF THE
      ! SGLS THAT DRAINED VIA CHANNEL-INCISION. SHOWS 1 FOR A LAKE THAT
      ! HAS STARTED TO DRAIN VIA CHANNEL INCISION AND 0 OTHERWISE. GETS
      ! MODIFIED IN THIS SUBROUTINE

      ! OUT_I,OUT_J: STORES THE COORDINATES OF CATCHMENT OUTLET. REMAINS
      ! UNMODIFIED IN THIS SUBROUTINE

      ! DEST_I,DEST_J: STORES COORDINATES OF THE DESTINATION OF THE
      ! MELTWATER PATH. REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! PATH_TOTAL: STORES TOTAL PATH LENGTH OF THE SUPRAGLACIAL
      ! MELTWATER CHANNEL. REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! PP,QQ: STORES COORDINATES OF LAKE NUCLEUS. REMAINS UNMODIFIED IN
      ! THIS SUBROUTINE

      ! K1,K2,K3,K4: CONTAINS EXTENT OF THE LAKE. REMAINS UNCHANGED IN
      ! THIS SUBROUTINE

      ! MCI: STORES THE LAKE. SHOWS 1 FOR LAKE PIXELS AND 0 OTHERWISE.
      ! REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! XIND,YIND: CONTAINS THE COORDINATES OF THE PIXELS LYING ALONG
      ! THE SUPRGLACIAL MELTWATER PATH. REMAINS UNMODIFIED IN THIS
      ! SUBROUTINE

      ! ROWS,COLS: ROWS AND COLUMNS OF THE STUDY AREA. REMAINS
      ! UNMODIFIED IN THIS SUBROUTINE

      ! SX: CONSTANT

      ! CMARK: SHOWS 1 FOR CREVASSED CELL AND 0 FOR NON CREVASSED CELL.
      ! UPDATED EVERY TIMESTEP IN THE CREV_CELLS1.F SUBROUTINE. REMAINS
      ! UNMODIFIED IN THIS SUBROUTINE

      ! SUMINC: STORES THE VOLUME OF WATER THAT HAS OVERFLOWN FROM A
      ! LAKE DURING CHANNEL INCISION. GETS MODIFIED IN THIS
      ! SUBROUTNE

      ! SURF: CONTAINS SURFACE HEIGHT (IN METERS) OF EVERY GRID CELL.
      ! REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! DEM: CONTAINS DEM SURFACE ELEVATION (IN METERS) FOR EVERY GRID
      ! CELL. REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! LAKE_HT: STORES THE REFERENCE LAKE HEIGHT. GETS CHANGED IN
      ! THIS SUBROUTINE

      ! LAKE_AR: STORES THE REFERENCE LAKE AREA. GETS CHANGED IN
      ! THIS SUBROUTINE

      ! LAKE_BASE: STORES THE REFERENCE LAKE VOLUME. GETS MODIFIED
      ! IN THIS SUBROUTINE

      ! CHN_DEC: STORES DEPTH OF SUPRAGLACIAL MELTWATER CHANNELS.
      ! GETS MODIFIED IN THIS SUBROUTINE

      ! RUNOFF_X: STORES THE RUNOFF THAT HAS BEEN TRANFERRED FROM A LAKE TO
      ! A NON-CREVASSED, NON-LAKE DEM CELL. THIS TRANSFER TAKES PLACE
      ! DURING LAKE OVERFLOW. GETS MODIFIED IN THIS SUBROUTINE

      ! RUNOFF_Y: STORES THE RUNOFF THAT HAS BEEN TRANFERRED FROM A LAKE TO
      ! A CREVASSED, NON-LAKE DEM CELL. THIS TRANSFER TAKES PLACE
      ! DURING LAKE OVERFLOW. GETS MODIFIED IN THIS SUBROUTINE

      ! RUNOFF_LAK: STORES THE RUNOFF THAT HAS BEEN TRANSFERRED FROM A LAKE
      ! TO ANOTHER EXISTING LAKE. THIS TRANSFER TAKES PLACE
      ! DURING LAKE OVERFLOW. GETS MODIFIED IN THIS SUBROUTINE

      ! RUNOFF: CONTAINS RUNOFF MELTWATER (IN METERS) FOR EVERY GRID
      ! CELL. REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! LAKE_RED: STORES THE HEIGHT OF WATER THAT HAS TO BE REDUCED FROM
      ! A LAKE AFTER IT HAS DRIANED VIA CHANNEL INCISION (IN EVERY
      ! TIMESTEP). GETS CHANGED IN THIS SUBROUTINE

      ! LAKE2RUN: STORES THE HEIGHT OF EXCESS LAKE WATER. THE EXCESS
      ! LAKE WATER IS A RESULT OF REDUCTION OF CHANNEL HEIGHT OF A LAKE
      ! DUE TO CHANNEL INCISION. GETS MODIFIED IN THIS SUBROUTINE

      ! LAK_NO_INC: STORES THE EXTENT OF THE LAKES THAT HAVE OVERFLOWN
      ! BUT WITHOUT CHANNEL INCISION. SHOWS 1 FOR LAKE PIXELS AND 0
      ! OTHERWISE. GETS MODIFIED IN THIS SUBROUTINE

      IMPLICIT NONE

      INTEGER :: SX,ROWS,COLS,MCI(COLS,ROWS),COUNT_C,A1,B1,A11,B11,PP,QQ
     &,T,I1,J1,C11,D11
     &,TIMESTEP,DRAIN_INCISE(COLS,ROWS),XIND(SX),YIND(SX),I,
     &CMARK(COLS,ROWS),NUMX,A,B,OUT_I,OUT_J,DEST_I,DEST_J,K1,K2,K3,K4

      DOUBLE PRECISION :: SLOPE,DEM(COLS,ROWS),PATH_TOTAL,
     &SURF(COLS,ROWS),LAKE_DEPTH(COLS,ROWS),VL,AR,
     &RK(COLS,ROWS),RUNOFF(COLS,ROWS),V_RUN,QIN,LAKE_HT(COLS,ROWS),
     &LAKE_AR(COLS,ROWS),LAKE_BASE(COLS,ROWS),PL,LAK_HT_REF,
     &LAK_AR_REF,LAK_HT_OLD,CHN_HT_OLD,GETA_OLD,CONST,ALPHA,BETA,CHN,
     &LAKE,GETA,V_N,GETA_NEW,DIFF,LAKE_RED(COLS,ROWS),CHN_DEC(COLS,ROWS)
     &,RUNOFF_X(COLS,ROWS),RUNOFF_Y(COLS,ROWS),RUNOFF_LAK(COLS,ROWS),
     &SUMINC(COLS,ROWS),LAKE2RUN(COLS,ROWS),CHN_HT_NEW,LAK_HT_NEW
     &,DIFF_CHN,V_NEW,SV,ELEV,SUMR

      LAKE_DEPTH=SURF-DEM
      SV=0

      VL=SUM(LAKE_DEPTH,MASK=MCI.GT.0)*100*100 ! LAKE VOLUME
      AR=SUM(MCI,MASK=MCI.GT.0)*100*100 ! LAKE AREA
      V_RUN=SUM(RUNOFF,MASK=MCI.GT.0)*100*100 ! TOTAL RUNOFF VOLUME
      COUNT_C=0

      IF (DRAIN_INCISE(PP,QQ).EQ.0) THEN
              LAKE_HT(PP,QQ)=SURF(PP,QQ)-DEM(PP,QQ)+RUNOFF(PP,QQ) ! REFERENCE HEIGHT
              LAKE_AR(PP,QQ)=AR ! REFERENCE LAKE AREA
              LAKE_BASE(PP,QQ)=VL+V_RUN ! REFERENCE LAKE VOLUME
      ELSE
      END IF

      NUMX=COUNT(XIND.NE.0)
      A1=OUT_I
      B1=OUT_J
      C11=XIND(NUMX-1)
      D11=YIND(NUMX-1)
      A11=DEST_I
      B11=DEST_J
      SLOPE=(DEM(A1,B1)-DEM(A11,B11))/PATH_TOTAL ! CHANNEL SLOPE

      QIN=0
      LAK_HT_REF=LAKE_HT(PP,QQ)
      LAK_AR_REF=LAKE_AR(PP,QQ)
      PL=LAK_AR_REF*LAK_HT_REF/LAKE_BASE(PP,QQ)

      LAK_HT_OLD=SURF(PP,QQ)-DEM(PP,QQ)+RUNOFF(PP,QQ) ! LAKE HEIGHT (PRESENT)
      CHN_HT_OLD=DEM(A1,B1)-DEM(PP,QQ) ! CHANNEL HEIGHT (PRESENT)
      GETA_OLD=LAK_HT_OLD-CHN_HT_OLD ! WATER LEVEL (PRESENT)

      CONST=((2*9.8)/(1+(0.25/4*SLOPE)))**1.5
      ALPHA=((0.25*1000)/(8*900.0*334000))*CONST
      BETA=CONST*((5*0.25)/(8*9.8*SLOPE))

      IF (CHN_HT_OLD.LT.0) CHN_HT_OLD=0
      IF (LAK_HT_OLD.LT.0) LAK_HT_OLD=0
      IF (GETA_OLD.LT.0) GETA_OLD=0

      ! ESTIMATING THE EVOLUTION OF CHANNEL AND LAKE DEPTH USING
      ! KINGSLAKE ET AL., 2015

      IF ((CHN_HT_OLD.GT.0).AND.(LAK_HT_OLD.GT.0)) THEN
              T=1
              CHN=CHN_HT_OLD
              LAKE=LAK_HT_OLD
              GETA=GETA_OLD
              V_N=VL+V_RUN ! LAKE VOLUME + LAKE RUNOFF

              DO WHILE ((LAKE.GT.0).AND.(CHN.GT.0).AND.(T.LE.86400)
     &)

              CHN_HT_NEW=CHN+(-ALPHA*(GETA**1.5))
              GETA_NEW=GETA+
     &1*((1/LAK_AR_REF)*((LAK_HT_REF/LAKE)**(PL-1))
     &*(QIN-BETA*GETA**1.5)+ALPHA*GETA**1.5)

              IF (GETA_NEW.LT.0) GETA_NEW=0
              LAK_HT_NEW=CHN_HT_NEW+GETA_NEW
              IF (CHN_HT_NEW.LT.0) CHN_HT_NEW=0
              IF (LAK_HT_NEW.LT.0) LAK_HT_NEW=0
              IF (V_NEW.LT.0) V_NEW=0

              CHN=CHN_HT_NEW
              LAKE=LAK_HT_NEW
              GETA=GETA_NEW
              V_N=V_NEW

              IF (T.LT.86400) T=T+1
              IF (T.EQ.86400) T=T+1
              END DO
              DIFF_CHN=CHN_HT_OLD-CHN_HT_NEW
              DIFF=LAK_HT_OLD-LAK_HT_NEW

              IF (GETA_NEW.GT.0) COUNT_C=1
      ELSE
              DIFF=0
              DIFF_CHN=0
      END IF

      IF (DIFF.GT.0) THEN
              IF (DRAIN_INCISE(PP,QQ).EQ.0) THEN ! FLAG THE LAKE WHEN IT OVERFLOWS
                      DRAIN_INCISE(PP,QQ)=1
              ELSE
              END IF

              ! CONVERTING EXTRA LAKE WATER INTO RUNOFF

              IF (GETA_NEW.GT.0) THEN
                      DO 140 I1=PP-K1,PP+K2
                      DO 150 J1=QQ-K3,QQ+K4
                      IF ((MCI(I1,J1).EQ.1)
     &.AND.(LAKE_DEPTH(I1,J1).GT.0)) THEN
                              IF (DIFF.GT.RUNOFF(PP,QQ)) THEN
                                      IF (SURF(I1,J1).GT.DEM(I1,J1)) 
     &THEN
                                              LAKE2RUN(I1,J1)=GETA_NEW
                                      END IF
                              ELSE
                                      IF (SURF(I1,J1).GT.DEM(I1,J1)) 
     &THEN
                                              LAKE2RUN(I1,J1)=DIFF_CHN
                                      END IF
                              END IF
                      END IF
150                   CONTINUE
140                   CONTINUE
              END IF

              ! STORING THE REDUCTION IN LAKE HEIGHT

              DO 701 I1=PP-K1,PP+K2
              DO 702 J1=QQ-K3,QQ+K4
              IF (MCI(I1,J1).EQ.1) THEN
                      LAKE_RED(I1,J1)=DIFF
              END IF
702           CONTINUE
701           CONTINUE

              ! STORING THE DECREASE IN CHANNEL ELEVATION

              NUMX=COUNT(XIND.GT.0)
              DO 420 I=1,NUMX-1
              A=XIND(I)
              B=YIND(I)
              IF (DIFF_CHN.GT.0) THEN
                      CHN_DEC(A,B)=CHN_DEC(A,B)+DIFF_CHN ! CHANNEL HEIGHT REDUCTION STORED
              END IF
420           CONTINUE

              ! ESTIMATING THE NEW VOLUME OF THE LAKE

              ELEV=DEM(PP,QQ)+LAK_HT_NEW
              DO 520 I1=PP-K1,PP+K2
              DO 530 J1=QQ-K3,QQ+K4
              IF (MCI(I1,J1).EQ.1) THEN
                      IF (ELEV.GE.DEM(I1,J1)) THEN
                              SV=SV+(ELEV-DEM(I1,J1)-CHN_DEC(I1,J1))
                      END IF
              END IF
530           CONTINUE
520           CONTINUE

              V_NEW=SV*10000
              SUMINC(PP,QQ)=(VL+V_RUN)-V_NEW ! WATER TRANSFERRED VIA
CHANNEL INCISION

              LAKE_DEPTH=SURF-DEM
              IF (LAKE_DEPTH(A11,B11).EQ.0) THEN
                      IF (CMARK(A11,B11).EQ.1) THEN
                              RUNOFF_Y(A11,B11)=RUNOFF_Y(A11,B11)+
     &SUMINC(PP,QQ)/10000
                      ELSE
                              RUNOFF_X(A11,B11)=RUNOFF_X(A11,B11)+
     &SUMINC(PP,QQ)/10000
                      END IF
              ELSE
                      RUNOFF_LAK(A11,B11)=RUNOFF_LAK(A11,B11)+
     &SUMINC(PP,QQ)/10000
              END IF
      ELSE
              SUMINC(PP,QQ)=0
      END IF

      RETURN

      END SUBROUTINE CHANNEL_INC1
