      SUBROUTINE FILL_DEPRESSION1(ROWS,COLS,DEM,RUNOFF,CMARK
     &,SURF,SX,SINK_I,SINK_J)

      ! THIS SUBROUTINE FILLS THE EMPTY SINK CELLS THEREBY FORMING LAKES

      ! DESCRIPTION OF VARIABLES:

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

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

      ! SINK_I,SINK_J: 1-D ARRAY. STORES THE NUCLEI OF THE LAKES. 
      ! REMAINS UNMODIFIED IN THIS SUBROUTINE

      ! SX: CONSTANT

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

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

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

      IMPLICIT NONE

      INTEGER :: SX,K1,K2,K3,K4,ROWS,COLS,CATCH_FLAG(COLS,ROWS),
     &MAX_LAK_EXT(COLS,ROWS),MARK(COLS,ROWS),XXX,I,J,OV_LAKE(COLS,ROWS)
     &,P,COUNTP,PP,QQ,XIND(SX),YIND(SX),XF(SX),YF(SX),I1,J1
     &,S,K,COUNTY,COUNTZ,A,B,N,II(9),JJ(9),ICLOW,NUM2,NUM1,KK,KKK,
     &SINK_I(SX),SINK_J(SX),VAL1(9),TEMP_XF,TEMP_YF,
     &MARK2(COLS,ROWS),MM,NN,INDI(2),III(9),TEMP_II,FLAG_SINK(COLS,ROWS)
     &,CMARK(COLS,ROWS),OUT_I,OUT_J,NTHREADS,OMP_GET_MAX_THREADS,
     &OMP_GET_THREAD_NUM,CHUNK,ISTART,IEND,TID,WAT_SINK(COLS,ROWS),NN2
     &,COND,TT,COUNTX,X,Y
      DOUBLE PRECISION :: DUM_LAK(COLS,ROWS),LAKE_DEPTH(COLS,ROWS),LOW,
     &SURF(COLS,ROWS),ACT_LAK(COLS,ROWS),RUN_O_LAK(COLS,ROWS),TEMP_SURF
     &,RUNOFF(COLS,ROWS),SURF_CATCH(SX),VALS(9),DEM(COLS,ROWS),
     &SURF_OUT,DEM_OUT,DIFF_H,MAX_LAKE_VOL,AV_RUNOFF_VOL,MAX_LAK_AR,
     &SUMR,R1,MAX_LAK_VOL,SURF1(COLS,ROWS),RUN1(COLS,ROWS),DEM_CATCH(SX)
     &,TEMP_DEM,HH,X_L,N_L(SX),ACT_LEV,FILL_VOL,NN1

      K1=0
      K2=0
      K3=0
      K4=0
      LAKE_DEPTH=SURF-DEM
      DUM_LAK=0.
      SURF_CATCH=0.
      VALS=0.
      MARK=0.
      CATCH_FLAG=0.
      MAX_LAK_EXT=0.
      XIND=0.
      YIND=0.
      XF=0.
      YF=0.
      VAL1=0.
      NN2=0
     

      XXX=COUNT(SINK_I.NE.0)

      ISTART=1
      IEND=XXX
     
!$OMP PARALLEL
!$OMP& DEFAULT(PRIVATE)
!$OMP& SHARED(SURF,DEM,RUNOFF,CMARK,SX,
!$OMP& XXX,NTHREADS,CHUNK,SINK_I,SINK_J,COLS,ROWS)
!$    NTHREADS=OMP_GET_MAX_THREADS()
!$    TID=OMP_GET_THREAD_NUM()
!$    CHUNK=(XXX+NTHREADS-1)/NTHREADS
!$    ISTART=TID*CHUNK+1
!$    IEND=MIN((TID+1)*CHUNK,XXX)

      DO 101 KKK=ISTART,IEND

      I=SINK_I(KKK)
      J=SINK_J(KKK)

      LAKE_DEPTH=SURF-DEM

      ! FINDING ALL THE NON LAKE SINK CELLS THAT HAVE NON-ZERO RUNOFF
      ! AND ARE EMPTY. THESE CELLS SHOULD HAVE POSITIVE SURFACE
      ! ELEVATION

      IF ((RUNOFF(I,J).GT.0.001).AND.(DEM(I,J).GT.0)
     &) THEN

      ! DETERMINING CORNER POINTS OF THE SEARCH BOX. THE SEARCH BOX WILL
      ! CONTAIN THE AREA WHERE WE SEARCH FOR ALL THE DEM CELLS THAT WILL
      ! POTENTIALLY CONTRIBUTE WATER INTO THE SINK CELL. THE SEARCH BOX
      ! IS A 200 BY 200 PIXEL BOX CENTERED AROUND THE SINK CELL

      PP=I
      QQ=J

      P=100
      COUNTP=1
      DO WHILE ((P.GE.0).AND.(COUNTP.EQ.1))
      IF ((PP-P).GE.2) THEN
              K1=P
              COUNTP=0
      ELSE
              P=P-1
              COUNTP=1
      END IF
      END DO

      P=100
      COUNTP=1
      DO WHILE ((P.GE.0).AND.(COUNTP.EQ.1))
      IF ((PP+P).LE.COLS-1) THEN
              K2=P
              COUNTP=0
      ELSE
              P=P-1
              COUNTP=1
      END IF
      END DO

      P=100
      COUNTP=1
      DO WHILE ((P.GE.0).AND.(COUNTP.EQ.1))
      IF ((QQ-P).GE.2) THEN
              K3=P
              COUNTP=0
      ELSE
              P=P-1
              COUNTP=1
      END IF
      END DO

      P=100
      COUNTP=1
      DO WHILE ((P.GE.0).AND.(COUNTP.EQ.1))
      IF ((QQ+P).LE.ROWS-1) THEN
              K4=P
              COUNTP=0
      ELSE
              P=P-1
              COUNTP=1
      END IF
      END DO

      XIND=XIND-XIND
      YIND=YIND-YIND
      SURF_CATCH=SURF_CATCH-SURF_CATCH
      XF=XF-XF
      YF=YF-YF
      CATCH_FLAG=0.
      
      ! DELINEATING THE SINK CELL'S CATCHMENT INSIDE THE SEARCH BOX

      DO 100 I1=PP-K1,PP+K2
      DO 110 J1=QQ-K3,QQ+K4
      S=I1
      K=J1
      NN=1
      COUNTY=1
      COUNTZ=1
      XIND=XIND-XIND
      YIND=YIND-YIND
      XF=XF-XF
      YF=YF-YF
      SURF_CATCH=SURF_CATCH-SURF_CATCH

      ! FINDING ALL THE DEM CELLS THAT WILL POTENTIALLY CONTRIBUTE
      ! RUNOFF INTO THE SINK CELL. THOSE CELLS THAT CAN POTENTIALLY
      ! CONTRIBUTE MELTWATER RUNOFF TO THE SINK CELL, WILL FORM A PART
      ! OF THE SINK CELL'S CATCHMENT

      DO WHILE ((NN.EQ.1).AND.(S.LE.COLS-1).AND.(K.LE.ROWS-1).AND.
     &(S.GE.2).AND.(K.GE.2))

      IF (COUNTZ.EQ.1) THEN
              XIND(COUNTY)=S
              YIND(COUNTY)=K
              COUNTY=COUNTY+1
              COUNTZ=COUNTZ+1
      END IF

      II=(/S-1,S-1,S-1,S,S,S,S+1,S+1,S+1/)
      JJ=(/K-1,K,K+1,K-1,K,K+1,K-1,K,K+1/)

      DO 120 N=1,9
      A=II(N)
      B=JJ(N)
      VALS(N)=DEM(A,B)
120   CONTINUE
      ICLOW=0
      LOW=999999

      DO 130 N=1,9
      IF (VALS(N).LT.LOW) THEN
              ICLOW=N
              LOW=VALS(N)
      END IF
130   CONTINUE

      IF (VALS(5).EQ.LOW) THEN
              NN=0
              S=II(5)
              K=JJ(5)
      ELSE
              S=II(ICLOW)
              K=JJ(ICLOW)
              NN=1
              XIND(COUNTY)=S
              YIND(COUNTY)=K
              COUNTY=COUNTY+1
      END IF
      END DO

      NUM2=COUNT(XIND.NE.0)

      ! ALL CELLS THAT HAVE THE SINK CELL AS THEIR RESPECTIVE POTENTIAL 
      ! DESTINATION ARE PART OF THE SINK CELL'S CATCHMENT

      IF ((PP.EQ.XIND(NUM2)).AND.(QQ.EQ.YIND(NUM2))) THEN
              NUM1=COUNT(XIND.NE.0)
              DO 140 KK=1,NUM1
              A=XIND(KK)
              B=YIND(KK)
              IF ((A.NE.0).AND.(B.NE.0)) THEN
                      CATCH_FLAG(A,B)=1 ! ACTIVATE THE CATCHMENT FLAG
              ELSE
                      CATCH_FLAG(A,B)=0
              END IF
140           CONTINUE
      ELSE
      END IF
110   CONTINUE
100   CONTINUE

      COUNTY=1

      ! IF A CATCHMENT EXISTS, THE NEXT STEP IS TO LOCATE THE OUTLET OF
      ! THE CATCHMENT. GENERALLY, THE OUTLET CELL IS DEFINED AS THE
      ! LOWEST LYING CATCHMENT'S BOUNDARY CELL

      ! LOCATING THE OUTLET CELL FOR THE DELINEATED CATCHMENT

      IF (SUM(CATCH_FLAG,MASK=CATCH_FLAG.EQ.1).GT.1) THEN ! IF LARGE CATCHMENT

              SURF_CATCH=SURF_CATCH-SURF_CATCH
              XF=XF-XF
              YF=YF-YF

              DO 150 I1=PP-K1,PP+K2
              DO 160 J1=QQ-K3,QQ+K4
              IF (CATCH_FLAG(I1,J1).EQ.1) THEN
                      II=(/I1-1,I1-1,I1-1,I1,I1,I1,I1+1,I1+1,I1+1/)
                      JJ=(/J1-1,J1,J1+1,J1-1,J1,J1+1,J1-1,J1,J1+1/)
                      DO 170 N=1,9
                      A=II(N)
                      B=JJ(N)
                      VAL1(N)=CATCH_FLAG(A,B)
170                   CONTINUE
                      IF (((VAL1(1).EQ.0).OR.(VAL1(2).EQ.0).OR.
     &(VAL1(3).EQ.0).OR.(VAL1(4).EQ.0).OR.(VAL1(6).EQ.0).OR.
     &(VAL1(7).EQ.0).OR.(VAL1(8).EQ.0).OR.(VAL1(9).EQ.0)).AND.
     &(VAL1(5).EQ.1)) THEN
                      A=II(5)
                      B=JJ(5)
                      XF(COUNTY)=A
                      YF(COUNTY)=B
                      SURF_CATCH(COUNTY)=DEM(A,B)
                      COUNTY=COUNTY+1 ! BOUNDARY CELLS OF CATCHMENT
                      END IF
              END IF
160           CONTINUE
150           CONTINUE

              ! SORTING THE ELEVATIONS IN INCREASING ORDER OF THEIR
              ! RESPECTIVE ELEVATIONS. BUBBLE SORT ALGORITHM IS USED FOR
              ! THAT PURPOSE

              NUM1=COUNT(SURF_CATCH.NE.0)
              DO 180 J1=0,NUM1-1
              DO 190 I1=1,NUM1-J1-1
              IF ((SURF_CATCH(I1).GT.SURF_CATCH(I1+1)).AND.
     &(SURF_CATCH(I1+1).GT.0)) THEN
              TEMP_SURF=SURF_CATCH(I1)
              SURF_CATCH(I1)=SURF_CATCH(I1+1)
              SURF_CATCH(I1+1)=TEMP_SURF
              TEMP_XF=XF(I1)
              XF(I1)=XF(I1+1)
              XF(I1+1)=TEMP_XF
              TEMP_YF=YF(I1)
              YF(I1)=YF(I1+1)
              YF(I1+1)=TEMP_YF
              END IF
190           CONTINUE
180           CONTINUE

              OUT_I=XF(1) ! COLUMN NO. OF OUTLET CELL
              OUT_J=YF(1) ! ROW NO. OF OUTLET CELL

              ! IN CASE THE OUTLET CELL IS A SINK CELL, A NEW OUTLET
              ! NEEDS TO BE LOCATED. THEREFORE, WE RUN THE NEW_OUTLET1
              ! SUBROUTINE

              IF ((OUT_I.GE.2).AND.(OUT_J.GE.2).AND.(OUT_I.LE.COLS-2)
     &.AND.(OUT_J.LE.ROWS-2)) THEN

              CALL NEW_OUTLET1(XF,YF,OUT_I,OUT_J,DEM,CATCH_FLAG,ROWS,
     &COLS,SX)

              END IF

              SURF_OUT=SURF(OUT_I,OUT_J) ! SURFACE ELEVATION OF OUTLET
              DEM_OUT=DEM(OUT_I,OUT_J) ! BASE ELEVATION OF OUTLET

              ! ISOLATING THE MAXIMUM EXTENT OF THE LAKE FROM OTHER
              ! LAKES

              DUM_LAK=0.
              MAX_LAK_EXT=0.

! WHEREVER THE CATCHMENT FLAG IS ENABLED, FIND THE MAXIMUM LAKE EXTENT
! AND CONSEQUENTLY THE MAXIMUM LAKE'S DEPTH DISTRIBUTION
              
              DO 403 I1=PP-K1,PP+K2
              DO 404 J1=QQ-K3,QQ+K4
              IF (CATCH_FLAG(I1,J1).EQ.1) THEN
                      DUM_LAK(I1,J1)=DEM_OUT-DEM(I1,J1)
                      IF (DUM_LAK(I1,J1).GE.0) THEN
                              MAX_LAK_EXT(I1,J1)=1
                      ELSE
                              DUM_LAK(I1,J1)=0
                              MAX_LAK_EXT(I1,J1)=0
                      END IF
              ELSE
                      DUM_LAK(I1,J1)=0
              END IF
404           CONTINUE
403           CONTINUE

              MAX_LAKE_VOL=10000*(SUM(DUM_LAK,MASK=DUM_LAK.GT.0)) ! MAXIMUM LAKE VOLUME
              AV_RUNOFF_VOL=10000*RUNOFF(PP,QQ) ! AVAILABLE RUNOFF VOLUME
              MAX_LAK_AR=10000*SUM(MAX_LAK_EXT,MASK=MAX_LAK_EXT.EQ.1) ! MAXIMUM LAKE AREA
              DIFF_H=(MAX_LAKE_VOL-AV_RUNOFF_VOL)/MAX_LAK_AR ! LAKE HEIGHT TO BE SUBTRACTED FROM MAXIMUM LAKE HEIGHT
             
              IF (SUM(MAX_LAK_EXT,MASK=CATCH_FLAG.EQ.1).EQ.0) THEN 
                      DIFF_H=0
              END IF

              IF (DIFF_H.GE.0) THEN ! CATCHMENT GETS HALF FULL OR JUST FULL

                      ! CATCHMENT PIXELS THAT MIGHT GET FILLED BY THE
                      ! LAKE ARE ISOLATED

                      COUNTY=0
                      SURF_CATCH=0.

                      DO 200 I1=PP-K1,PP+K2
                      DO 210 J1=QQ-K3,QQ+K4                      
                      IF ((MAX_LAK_EXT(I1,J1).EQ.1)
     &.AND.(SURF(I1,J1).LE.DEM_OUT).AND.(SURF(I1,J1).GT.0)) THEN
                              COUNTY=COUNTY+1
                              SURF_CATCH(COUNTY)=SURF(I1,J1)
                      END IF
210                   CONTINUE
200                   CONTINUE

                      ! ISOLATED CATCHMENT ELEVATIONS ARRANGED IN THE
                      ! INCREASING ORDER

                      NUM1=COUNT(SURF_CATCH.NE.0)

                      DO 413 J1=0,NUM1-1
                      DO 414 I1=1,NUM1-J1-1
                      IF ((SURF_CATCH(I1).GT.SURF_CATCH(I1+1))
     &.AND.(SURF_CATCH(I1+1).GT.0)) THEN
                              TEMP_SURF=SURF_CATCH(I1)
                              SURF_CATCH(I1)=SURF_CATCH(I1+1)
                              SURF_CATCH(I1+1)=TEMP_SURF
                      END IF
414                   CONTINUE
413                   CONTINUE

                      ! FIND THE WATER SURFACE ELEVATION USING AN
                      ! ITERATIVE PROCESS

                      HH=SURF(PP,QQ) ! INITIAL GUESS OF WATER LEVEL
                      COND=1
                      DO WHILE ((COND.EQ.1).AND.(HH.LT.DEM_OUT))
                      
                      ! ESTIMATING WATER DEPTHS IN THE DEPRESSION FOR HH

                      N_L=0.
                      X_L=0
                      DO 405 I1=1,NUM1
                      IF ((HH-SURF_CATCH(I1)).GT.0) THEN
                              N_L(I1)=HH-SURF_CATCH(I1)
                      END IF
405                   CONTINUE

                      X_L=SUM(N_L) ! TOTAL LAKE DEPTHS
                      FILL_VOL=X_L*10000 ! VOLUME OF FILLING
                      
                      ! IF FILLING VOLUME LESS THAN WHAT IS AVAILABLE,
                      ! INCREASE THE WATER SURFACE ELEVATION AND
                      ! CONTINUE THE PROCESS

                      IF (FILL_VOL.LT.AV_RUNOFF_VOL) THEN
                              IF ((RUNOFF(PP,QQ).GT.0.001).AND.
     &(RUNOFF(PP,QQ).LE.0.01)) THEN
                                      HH=HH+0.0001
                              END IF

                              IF ((RUNOFF(PP,QQ).GT.0.01).AND.
     &(RUNOFF(PP,QQ).LE.0.1)) THEN
                                      HH=HH+0.001
                              END IF

                              IF ((RUNOFF(PP,QQ).GT.0.1).AND.
     &(RUNOFF(PP,QQ).LE.1)) THEN
                                      HH=HH+0.01
                              END IF
                              IF ((RUNOFF(PP,QQ).GT.1)) THEN
                                      HH=HH+0.1
                              END IF
                              
                              COND=1   
                      ELSE
                              COND=0
                      END IF
                      END DO

                      ! NEW WATER SURFACE ELEVATION FOUND

                              IF (HH.GT.DEM_OUT) THEN
                                     HH=DEM_OUT
                              END IF
                              ACT_LEV=HH

                      ! ALL THOSE CATCHMENT CELLS WITH SURFACE ELEVATION
                      ! LOWER THAN OR EQUAL TO 'ACT_LEV' GET FILLED I.E.
                      ! THEIR RESPECTIVE SURFACE ELEVATIONS GET RAISED
                      ! TO ACT_LEV

!$OMP CRITICAL                      
                      DO 406 I1=PP-K1,PP+K2
                      DO 407 J1=QQ-K3,QQ+K4
                      IF ((SURF(I1,J1).LE.ACT_LEV).AND.
     &(MAX_LAK_EXT(I1,J1).EQ.1)) THEN
                              SURF(I1,J1)=ACT_LEV
                      END IF
407                   CONTINUE
406                   CONTINUE

                      RUNOFF(PP,QQ)=0
!$OMP END CRITICAL
                      CATCH_FLAG=0.
                      ACT_LEV=0
                      HH=0
                      N_L=0.

              ELSE ! CATCHMENT READY TO OVERFLOW

                      ! AS THE LAKE IS OVERFLOWING, ACTIVATE THE OV_LAKE
                      ! MATRIX FOR THE LAKE
!$OMP CRITICAL
                      DO 250 I1=PP-K1,PP+K2
                      DO 260 J1=QQ-K3,QQ+K4                  
                      IF (MAX_LAK_EXT(I1,J1).EQ.1) THEN
                              SURF(I1,J1)=DEM_OUT
                      END IF
260                   CONTINUE
250                   CONTINUE

                      RUNOFF(PP,QQ)=0
                      RUNOFF(PP,QQ)=(ABS(AV_RUNOFF_VOL-MAX_LAKE_VOL))
     &/10000 ! REST WILL OVERFLOW

!$OMP END CRITICAL                      
                      
              END IF
              MAX_LAK_EXT=0. ! INITIALISING THE MAXIMUM LAKE EXTENT MATRIX
              CATCH_FLAG=0. ! CLEARING THE CATCHMENT
              DUM_LAK=0. ! CLEARING MAXIMUM LAKE DEPTH DISTRIBUTION
              DIFF_H=0.0 ! CLEARING THE SCALAR VARIABLE
              MARK=0. ! CLEARING THE LAKE DELINEATION MATRIX
              SURF_CATCH=0.
      END IF

      IF (SUM(CATCH_FLAG).EQ.1) THEN ! IF SMALL CATCHMENT!
              INDI=MAXLOC(CATCH_FLAG,MASK=CATCH_FLAG.GT.0) ! DELINEATE
CATCHMENT
              MM=INDI(1)
              NN=INDI(2)
              II=(/MM-1,MM-1,MM-1,MM,MM,MM,MM+1,MM+1,MM+1/)
              JJ=(/NN-1,NN,NN+1,NN-1,NN,NN+1,NN-1,NN,NN+1/)
              DO 220 N=1,9
              VALS(N)=DEM(II(N),JJ(N))
              III(N)=N
220           CONTINUE

              DO 230 J1=0,9-1
              DO 240 I1=1,9-J1-1
              IF (VALS(I1).GT.VALS(I1+1)) THEN
                      TEMP_SURF=VALS(I1)
                      VALS(I1)=VALS(I1+1)
                      VALS(I1+1)=TEMP_SURF
                      TEMP_II=III(I1)
                      III(I1)=III(I1+1)
                      III(I1+1)=TEMP_II
              END IF
240           CONTINUE
230           CONTINUE
              OUT_I=II(III(2))
              OUT_J=JJ(III(2))
              DUM_LAK=0.
              AV_RUNOFF_VOL=0.0
              MAX_LAKE_VOL=0.
              DIFF_H=0.0
!$OMP CRITICAL
              DUM_LAK(PP,QQ)=SURF(OUT_I,OUT_J)-DEM(PP,QQ) ! MAX LAKE DEPTH
              AV_RUNOFF_VOL=RUNOFF(PP,QQ)*100*100 ! AVAILABLE RUNOFF
              MAX_LAKE_VOL=DUM_LAK(PP,QQ)*100*100 ! MAX LAKE VOLUME
              DIFF_H=(MAX_LAKE_VOL-AV_RUNOFF_VOL)/10000 ! LAKE HEIGHT TO BE SUBTRACTED FROM MAXIMUM LAKE HEIGHT
              IF (DIFF_H.GE.0) THEN

                      SURF(PP,QQ)=DEM(OUT_I,OUT_J)-DIFF_H
                      RUNOFF(PP,QQ)=0
              ELSE
                      SURF(PP,QQ)=SURF(OUT_I,OUT_J) 
                      RUNOFF(PP,QQ)=(ABS(AV_RUNOFF_VOL-MAX_LAKE_VOL))
     &/10000
              END IF          
              MAX_LAK_EXT=0.
              CATCH_FLAG=0.
              DUM_LAK=0.
              DIFF_H=0.
!$OMP END CRITICAL              
      END IF

      END IF
101   CONTINUE
!$OMP END PARALLEL
      RETURN

      END SUBROUTINE FILL_DEPRESSION1

