SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,
     1     NDOFEL,NRHS,NSVARS,PROPS,NPROPS,COORDS,
     2     MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
     3     KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,
     4     ADLMAG,PREDEF,NPREDF,LFLAGS,MLVARX,
     5     DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,
     6     PERIOD)
C     
      INCLUDE 'ABA_PARAM.INC'
      PARAMETER ( ZERO = 0.D0, HALF = 0.5D0, 
     1     ONE = 1.D0, SEVEN=7.0D0, EIGHT=8.0D0 )
C
      DIMENSION RHS(MLVARX,*),
     1     AMATRX(NDOFEL,NDOFEL),
     2     SVARS(NSVARS),ENERGY(8),PROPS(*),
     3     COORDS(MCRD,NNODE),U(NDOFEL),
     4     DU(MLVARX,*),V(NDOFEL),A(NDOFEL),
     5     TIME(2),PARAMS(3),JDLTYP(MDLOAD,*),
     6     ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     7     PREDEF(2,NPREDF,NNODE),LFLAGS(*),
     8     JPROPS(*)
C
C     USER ELEMENT TO MODEL THE HYDRODYNAMIC 
C     PRESSURE OF RESERVOIR ON DAM USING 
C     WESTERGAARD ADDED MASS TECHNIQUE. 
C     
C     PROPERTIES:
C           PROPS(1) = Depth of reservoir
C           PROPS(2) = Y coordinate of water level
C           PROPS(3) = Density of water
C
      RESERVOIR_DEPTH    = PROPS(1)
      YCOORD_WATER = PROPS(2)
      WATER_DENSITY  = PROPS(3)
      YCOORD_ELCENTROID = 
     *    HALF * (COORDS(2,2) + COORDS(2,1))
      ELCENTROID_DEPTH = 
     *    YCOORD_WATER - YCOORD_ELCENTROID
      RHO = (SEVEN/EIGHT) * WATER_DENSITY *
     *    SQRT(RESERVOIR_DEPTH * ELCENTROID_DEPTH)
C     
      ALEN = ABS(COORDS(2,2)-COORDS(2,1))
      AM   = HALF*RHO*ALEN
C
      DO K1 = 1, NDOFEL                      
        DO KRHS = 1, NRHS
          RHS(K1,KRHS) = ZERO
        END DO
        DO K2 = 1, NDOFEL
          AMATRX(K2,K1) = ZERO
        END DO
      END DO

      do k1 = 1, nsvars
         svars(k1) = zero
      end do
C
      IF (LFLAGS(3).EQ.1) THEN
C       Normal incrementation
        IF (LFLAGS(1).EQ.11
     *        .OR.LFLAGS(1).EQ.12) THEN
C        *DYNAMIC
          BETA  = PARAMS(2)
          DADU = ONE/(BETA*DTIME**2)
C     
          AMATRX(1,1) = AM*DADU
          AMATRX(3,3) = AM*DADU
          RHS(1,1) = RHS(1,1)-AM*A(1)
          RHS(3,1) = RHS(3,1)-AM*A(3)
          ENERGY(1)= HALF*AM*(V(1)*V(1)+V(3)*V(3))
        END IF
      ELSE IF (LFLAGS(3).EQ.4) THEN
C       Mass matrix
         AMATRX(1,1) = AM
         AMATRX(3,3) = AM
      ELSE IF (LFLAGS(3).EQ.5) THEN
C       Half-step residual calculation
         RHS(1,1) = RHS(1,1)-AM*A(1)
         RHS(3,1) = RHS(3,1)-AM*A(3)
      ELSE IF (LFLAGS(3).EQ.6) THEN
C       Initial acceleration calculation
         AMATRX(1,1) = AM
         AMATRX(3,3) = AM
         ENERGY(1) = HALF*AM*(V(1)*V(1)+V(3)*V(3))
      END IF
C
      RETURN
      END