      SUBROUTINE STRESS_INT1(CREV_DEP,RUNOFF,WAT_DEP,STRESS,ICE,C_MARK
     &,I,J,COLS,ROWS)

      ! THIS SUBROUITNE CALCULATES THE STRESS INTENSITY FACTOR FOR EVERY
      ! WATER-FILLED NON-LAKE CREVASSE

      ! DESCRIPTION OF VARIABLES:

      ! 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

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

      ! CREV_DEP: CONTAINS CREVASSE DEPTH FOR EVERY GRID CELL.
      ! GETS MODIFIED IN THIS SUBROUTINE

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

      ! STRESS: CONTAINS STRESS (IN kPa) FOR EVERY GRID CELL. REMAINS
      ! UNMODIFIED IN THIS SUBROUTINE

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

      ! WAT_DEP: SHOWS DEPTH OF WATER (IN METERS) IN CREVASSED CELLS. FOR
      ! NON-CREVASSED CELLS, THE VALUE IS 0. GETS MODIFIED IN THIS
      ! SUBROUTINE

      ! I,J: COORDINATES OF A GRID CELL. UNMODIFIED IN THIS SUBROUTINE

      IMPLICIT NONE
      INTEGER :: COLS,ROWS,TIMESTEP,C_MARK(COLS,ROWS),I,J
      DOUBLE PRECISION :: CREV_DEP(COLS,ROWS),RUNOFF(COLS,ROWS),
     &WAT_DEP(COLS,ROWS),STRESS(COLS,ROWS),ICE(COLS,ROWS),KI,
     &WID,V_EMPT,REM_HT,V_RUN,ROOT_FIN,KIC

      WID=0.6
      KIC=150000
      ! IF RUNOFF>0, ADD WATER, CALCULATE KI AND UPDATE CREVASSE DEPTH
      ! ITERATIVELY
      IF (RUNOFF(I,J).GT.0) THEN
              ! ADD WATER TILL ALL MELTWATER IS ACCOMMODATED IN THE
              ! CREVASSE
              DO WHILE ((RUNOFF(I,J).GT.0)
     &.AND.(CREV_DEP(I,J).LT.ICE(I,J)))

              REM_HT= CREV_DEP(I,J)-WAT_DEP(I,J) ! PART OF CREVASSE NOT FILLED
              V_EMPT=REM_HT*WID*100   ! VOLUME OF EMPTY SPACE
              V_RUN=RUNOFF(I,J)*100*100    ! RUNOFF AVAILABLE IN THE CELL
              IF (V_RUN.GT.V_EMPT) THEN ! RUNOFF>EMPTY SPACE VOLUME
                      WAT_DEP(I,J)=CREV_DEP(I,J)  ! CREVASSE FILLED
                      RUNOFF(I,J)=RUNOFF(I,J)-V_EMPT/(100*100) ! RUNOFF UPDATED
              ELSE IF (V_RUN.EQ.V_EMPT) THEN ! RUNOFF=EMPTY SPACE VOLUME
                      WAT_DEP(I,J)=CREV_DEP(I,J) ! CREVASSE FILLED
                      RUNOFF(I,J)=0         ! RUNOFF EMPTIED
              ELSE IF (V_RUN.LT.V_EMPT) THEN ! RUNOFF<EMPTY SPACE VOLUME
                      WAT_DEP(I,J)=WAT_DEP(I,J)+V_RUN/(WID*100) ! WATER DEPTH INC.
                      RUNOFF(I,J)=0
              ELSE
              END IF
              ! CALCULATE KI
              KI=(1.12*1000*STRESS(I,J)*SQRT(3.141*CREV_DEP(I,J)))
     &-(0.683*900*9.8*
     &CREV_DEP(I,J)**1.5)+(0.683*1000*9.8*(WAT_DEP(I,J))**1.5)
              ! IF KI>KIC ITERATIVELY INCREASE THE CREVASSE DEPTH
              IF (KI.GT.KIC) THEN
                      DO WHILE ((KI.GT.KIC)
     &.AND.(CREV_DEP(I,J).LT.ICE(I,J)))
                      CREV_DEP(I,J)=CREV_DEP(I,J)+0.1 ! UPDATE CREVASSE DEPTH
                      ! RECALCULATE KI
                      KI=(1.12*1000*STRESS(I,J)
     &*SQRT(3.141*CREV_DEP(I,J)))
     &-(0.683*900*9.8*CREV_DEP(I,J)**1.5)
     &+(0.683*1000*9.8*(WAT_DEP(I,J))**1.5)
                      END DO
                      IF (CREV_DEP(I,J).GE.ICE(I,J)) 
     &CREV_DEP(I,J)=ICE(I,J) ! MOULIN FORMED
              END IF

              IF ((CREV_DEP(I,J).EQ.WAT_DEP(I,J))
     &.AND.(RUNOFF(I,J).GT.0)) THEN
                      CREV_DEP(I,J)=CREV_DEP(I,J)+0.01
                      KI=(1.12*1000*STRESS(I,J)
     &*SQRT(3.141*CREV_DEP(I,J)))
     &-(0.683*900*9.8*CREV_DEP(I,J)**1.5)
     &+(0.683*1000*9.8*(WAT_DEP(I,J))**1.5)
              END IF
              END DO
      ELSE
              ! IF NO RUNOFF CALCULATE KI, UPDATE CREVASSE DEPTH TILL
              ! KI>KIC
              KI=(1.12*1000*STRESS(I,J)*
     &SQRT(3.141*CREV_DEP(I,J)))-(0.683*900*9.8*
     &CREV_DEP(I,J)**1.5)+(0.683*1000*9.8*(WAT_DEP(I,J))**1.5)
              IF (KI.GT.KIC) THEN
                      DO WHILE ((KI.GT.KIC)
     &.AND.(CREV_DEP(I,J).LT.ICE(I,J)))
                      CREV_DEP(I,J)=CREV_DEP(I,J)+0.1
                      ! RECALCULATE KI
                      KI=(1.12*1000
     &*STRESS(I,J)*SQRT(3.141*CREV_DEP(I,J)))
     &-(0.683*900*9.8*CREV_DEP(I,J)**1.5)
     &+(0.683*1000*9.8*(WAT_DEP(I,J))**1.5)
                      END DO
                      IF (CREV_DEP(I,J).GE.ICE(I,J)) 
     &CREV_DEP(I,J)=ICE(I,J)
              END IF
      END IF

      RETURN
      END SUBROUTINE STRESS_INT1
