SUBROUTINE FRIC(LM,TAU,DDTDDG,DDTDDP,DSLIP,SED,SFD,
     1     DDTDDT,PNEWDT,STATEV,DGAM,TAULM,PRESS,DPRESS,DDPDDH,
     2     SLIP,KSTEP,KINC,TIME,DTIME,NOEL,CINAME,SLNAME,
     3     MSNAME,NPT,NODE,NPATCH,COORDS,RCOORD,DROT,TEMP,
     4     PREDEF,NFDIR,MCRD,NPRED,NSTATV,CHRLNGTH,PROPS,NPROPS)
C                    
      INCLUDE 'ABA_PARAM.INC'
C                    
      CHARACTER*80 CINAME,SLNAME,MSNAME
      DIMENSION TAU(NFDIR),DDTDDG(NFDIR,NFDIR),DDTDDP(NFDIR),
     1     DSLIP(NFDIR),DDTDDT(NFDIR,2),STATEV(*),
     2     DGAM(NFDIR),TAULM(NFDIR),SLIP(NFDIR),
     3     COORDS(MCRD),RCOORD(MCRD),DROT(2,2),TEMP(2),
     4     PREDEF(2,*),TIME(2),PROPS(NPROPS)
C                    
      DIMENSION GAMMA(2),TAURV(2),GAMRV(2)
C                    
      PARAMETER(ZERO=0.0D0)
      REAL FINC,PM,DSLIPEQ,DGAMEQ,CPT
C	  
C     IMPLEMENTATION OF 3-D COULOMB FRICTION USING PENALTY METHOD
C     
C     VARIABLES USED:
C        XMU = COEFFICIENT OF FRICTION
C        GCRIT = CRITICAL ELASTIC SLIP
C        STIFF = ARTIFICIAL STIFFNESS
C        GAMMA(1:2) = TOTAL SLIP
C        STATEV(1:2) = ELASTIC SLIP
C        TAUCRIT = CRITICAL FRICTIONAL STRESS*
C
C      NODEID(1) = 304
C	  
	  A  = PROPS(1)
	  B  = PROPS(2)
	  C  = PROPS(3)
	  D  = PROPS(4)
	  E  = PROPS(5)
C                      
      IF (LM .EQ. 2) THEN 
  	      RETURN
      ENDIF	
C 	  
      IF ((KSTEP .EQ. 1) .AND. (KINC .EQ. 1)) THEN 
	      CPT = ZERO
      ENDIF
C	  
C	  CPT = 1
C	  
      IF (CPT .EQ. 1) THEN
 	     CPT = 1
      ENDIF
C	  
      LM = 0
C	  
      IF (PRESS .LT. ZERO) THEN 
          RETURN
C           STATEV(1)   = ZERO
C           STATEV(2)   = ZERO
C           GAMMA(1)    = ZERO
C           GAMMA(2)    = ZERO
C 		  DDTDDG(1,1) = 1000000
C 		  DDTDDG(1,2) = ZERO
C 		  DDTDDG(2,1) = ZERO
C 		  DDTDDG(2,2) = 1000000
C 		  DDTDDP(1)   = ZERO
C 		  DDTDDP(2)   = ZERO
C           TAU(1)      = ZERO
C           TAU(2)      = ZERO
C           DSLIP(1)    = ZERO
C           DSLIP(2)    = ZERO  
      ELSE
C	  
	  GAMMA(1) = STATEV(1) + DGAM(1)
	  GAMMA(2) = STATEV(2) + DGAM(2) 		      
C
      DGAMEQ = (DGAM(1) + DGAM(2))/2
	  DSLIPEQ = (DSLIP(1)+ DSLIP(2))/2
C	  
	      IF (DGAMEQ*DSLIPEQ .GT. ZERO) THEN  	   
			  CPT = 1
			  PM = PRESS
			  TAURV(1) = TAU(1)
			  TAURV(2) = TAU(2)
			  GAMRV(1) = GAMMA(1)
			  GAMRV(2) = GAMMA(2)
			  
			  DDTDDG(1,1) = ZERO
			  DDTDDG(1,2) = ZERO	
			  DDTDDG(2,1) = ZERO 		 
			  DDTDDG(2,2) = ZERO
			  DDTDDP(1)   = TAURV(1)
			  DDTDDP(2)   = TAURV(2)
C				  
			  TAU(1) = TAURV(1)*PRESS
			  TAU(2) = TAURV(2)*PRESS
C			  
			  DSLIP(1)  = DGAM(1)
			  DSLIP(2)  = DGAM(2)
C	C	C	C	  
			  STATEV(1) = GAMMA(1)
			  STATEV(2) = GAMMA(2)
C              RETURN
C		  
          ELSE
		      IF (CPT .EQ. ZERO) THEN
C	C	C	C		   
				  DDTDDG(1,1) = A*B*EXP(-B*GAMMA(1))*PRESS
				  DDTDDG(1,2) = ZERO	
				  DDTDDG(2,1) = ZERO 		 
				  DDTDDG(2,2) = A*B*EXP(-B*GAMMA(2))*PRESS
				  DDTDDP(1)   = A*(1-EXP(-B*GAMMA(1)))
				  DDTDDP(2)   = A*(1-EXP(-B*GAMMA(2)))
C				  
				  TAU(1) = (A*(1-EXP(-B*GAMMA(1))))*PRESS
				  TAU(2) = (A*(1-EXP(-B*GAMMA(2))))*PRESS
C	C	C	C
				  DSLIP(1)  = DGAM(1)
				  DSLIP(2)  = DGAM(2)
C	C	C	C	  
				  STATEV(1) = GAMMA(1)
				  STATEV(2) = GAMMA(2)
C	C		  
              ELSEIF (CPT .EQ. 1) THEN
C	C		  	   
				  DDTDDG(1,1) = C*PRESS
				  DDTDDG(1,2) = ZERO	
				  DDTDDG(2,1) = ZERO 		 
				  DDTDDG(2,2) = C*PRESS
				  DDTDDP(1)   = (C*(GAMMA(1)-GAMRV(1))+TAURV(1))
				  DDTDDP(2)   = (C*(GAMMA(2)-GAMRV(2))+TAURV(2))
C				  
				  TAU(1) = (C*(GAMMA(1)-GAMRV(1))+TAURV(1))*PRESS
				  TAU(2) = (C*(GAMMA(2)-GAMRV(2))+TAURV(2))*PRESS
C	C	C	C
				  DSLIP(1)  = DGAM(1)
				  DSLIP(2)  = DGAM(2)
C	C	C	C	  
				  STATEV(1) = GAMMA(1)
				  STATEV(2) = GAMMA(2)
C              ELSE
              ENDIF			  
          ENDIF
      ENDIF
      RETURN
      END