!DEC$ FREEFORM

    SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,&
    DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,REDEF,DRED,&
    CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,&
    CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
    
! ----------------------------------------------------------------
!    umat for compressible neo-hookean hyperelasticity
!    cannot be used for plane stress
! ----------------------------------------------------------------

    USE ABQinterface
    
    IMPLICIT NONE
    
    !---------------------variables from Abaqus interface----------------------
    CHARACTER(LEN=80):: CMNAME
    INTEGER(KIND=AbqIK)::NOEL,NPT,LAYER,KSPT,KINC,NDI,NSHR,NTENS,NSTATV,NPROPS
    INTEGER(KIND=AbqIK),DIMENSION(4):: JSTEP
    REAL(KIND=AbqRK):: DTIME,TEMP,DTEMP,SSE,SPD,SCD,CELENT,DRPLDT,PNEWDT,RPL
    REAL(KIND=AbqRK),DIMENSION(NTENS):: STRESS,STRAN,DSTRAN
    REAL(KIND=AbqRK),DIMENSION(NSTATV):: STATEV
    REAL(KIND=AbqRK),DIMENSION(NTENS,NTENS):: DDSDDE
    REAL(KIND=AbqRK),DIMENSION(NTENS)::DDSDDT
    REAL(KIND=AbqRK),DIMENSION(NTENS):: DRPLDE
    REAL(KIND=AbqRK),DIMENSION(2):: TIME
    REAL(KIND=AbqRK),DIMENSION(1):: REDEF,DRED
    REAL(KIND=AbqRK),DIMENSION(NPROPS):: PROPS
    REAL(KIND=AbqRK),DIMENSION(3):: COORDS
    REAL(KIND=AbqRK),DIMENSION(3,3):: DROT,DFGRD0,DFGRD1

!-------------------------------------local arrays---------------------------
!    B          - deviatoric right cauchy-green tensor
!    BP         - principal values of B
!    BN         - principal direction of B
!    distortion - deviatoric deformation gradient (distortion tensor)
! ----------------------------------------------------------------------------
     REAL(KIND=AbqRK):: B(6), BP(3), BN(3, 3), distortion(3,3)
     REAL(KIND=AbqRK):: E, NU, E2, E3, R, E4, C10, D1, determinant, scaling, TB
     INTEGER(KIND=AbqIK):: i, j
     
!    PROPS(1) - E
!    PROPS(2) - NU
! ----------------------------------------------------------------------------
!    elastic Properties
! ----------------------------------------------------------------------------
      E=PROPS(1)      
      NU=PROPS(2)      
      C10=E/(4.0_AbqRK*(1.0_AbqRK+NU))      
      D1=6.0_AbqRK*(1.0_AbqRK-2.0_AbqRK*NU)/E

!    jacobian and distortion tensor

      determinant=DFGRD1(1, 1)*DFGRD1(2, 2)*DFGRD1(3, 3)&
        -DFGRD1(1, 2)*DFGRD1(2, 1)*DFGRD1(3, 3)
      IF(NSHR.EQ.3) THEN        
         determinant=determinant+DFGRD1(1, 2)*DFGRD1(2, 3)*DFGRD1(3, 1)&     
              +DFGRD1(1, 3)*DFGRD1(3, 2)*DFGRD1(2, 1)&     
              -DFGRD1(1, 3)*DFGRD1(3,1)*DFGRD1(2, 2)&     
              -DFGRD1(2, 3)*DFGRD1(3, 2)*DFGRD1(1, 1)   
      END IF      
      scaling=determinant**(-1.0_AbqRK/3.0_AbqRK)      
      DO i=1, 3        
        DO j=1, 3          
          distortion(j, i)=scaling*DFGRD1(j, i)        
        END DO
      END DO
!    calculate deviatoric left cauchy-green deformation tensor
      
      B(1)=distortion(1, 1)**2+distortion(1, 2)**2+distortion(1, 3)**2      
      B(2)=distortion(2, 1)**2+distortion(2, 2)**2+distortion(2, 3)**2      
      B(3)=distortion(3, 3)**2+distortion(3, 1)**2+distortion(3, 2)**2      
      B(4)=distortion(1, 1)*distortion(2, 1)+distortion(1, 2)*distortion(2, 2)&     
             +distortion(1, 3)*distortion(2, 3)      
      IF(NSHR.EQ.3) THEN        
        B(5)=distortion(1, 1)*distortion(3, 1)+distortion(1, 2)*distortion(3, 2)&     
              +distortion(1, 3)*distortion(3, 3)        
        B(6)=distortion(2, 1)*distortion(3, 1)+distortion(2, 2)*distortion(3, 2)&     
              +distortion(2, 3)*distortion(3, 3)      
      END IF

!    calculate the stress

      TB=(B(1)+B(2)+B(3))/3.0_AbqRK      
      E2=2.0_AbqRK*C10/determinant      
      E3=2.0_AbqRK/D1*(2.0_AbqRK*determinant-1.0_AbqRK)      
      R=2.0_AbqRK/D1*(determinant-1.0_AbqRK)      
      DO i=1,NDI      
        STRESS(i)=E2*(B(i)-TB)+R      
      END DO
      DO i=NDI+1,NDI+NSHR        
        STRESS(i)=E2*B(i)      
      END DO
!    calculate the stiffness

      E4=E2*2.0_AbqRK/3.0_AbqRK
          
      DDSDDE(1, 1)= E4*(B(1)+TB)+E3      
      DDSDDE(2, 2)= E4*(B(2)+TB)+E3      
      DDSDDE(3, 3)= E4*(B(3)+TB)+E3      
      DDSDDE(1, 2)=-E4*(B(1)+B(2)-TB)+E3      
      DDSDDE(1, 3)=-E4*(B(1)+B(3)-TB)+E3      
      DDSDDE(2, 3)=-E4*(B(2)+B(3)-TB)+E3      
      DDSDDE(1, 4)= E4*B(4)/2.0_AbqRK      
      DDSDDE(2, 4)= E4*B(4)/2.0_AbqRK      
      DDSDDE(3, 4)=-E4*B(4)      
      DDSDDE(4, 4)= E2*(B(1)+B(2))/2.0_AbqRK
      IF(NSHR.EQ.3) THEN      
        DDSDDE(1, 5)= E4*B(5)/2.0_AbqRK        
        DDSDDE(2, 5)=-E4*B(5)        
        DDSDDE(3, 5)= E4*B(5)/2.0_AbqRK        
        DDSDDE(1, 6)=-E4*B(6)        
        DDSDDE(2, 6)= E4*B(6)/2.0_AbqRK        
        DDSDDE(3, 6)= E4*B(6)/2.0_AbqRK        
        DDSDDE(5, 5)= E2*(B(1)+B(3))/2.0_AbqRK
        DDSDDE(6, 6)= E2*(B(2)+B(3))/2.0_AbqRK        
        DDSDDE(4,5)= E2*B(6)/2.0_AbqRK        
        DDSDDE(4,6)= E2*B(5)/2.0_AbqRK        
        DDSDDE(5,6)= E2*B(4)/2.0_AbqRK      
      END IF      
      DO i=1, NTENS      
        DO j=1, i-1        
          DDSDDE(i, j)=DDSDDE(j, i)        
        END DO     
      END DO
      
      PNEWDT=1.5 !set PNEWDT to a value above 1.0, to indicate that the time
                 !increment can be increased
      
      RETURN      
      END SUBROUTINE UMAT
