cohesive model of Needleman was input into abaqus, but wrong???

 I have written  vumat codes of a cohesive model of Needleman (whose paper entilted "Numerical simulations of fast crack growth in brittle solids"), and I have implemented it into ABAQUS/Explicit, but the predicted result is wrong. I don't know why. The corresponding fortran file is as follows:

      subroutine vumat(

C Read only (unmodifiable)variables -

     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,

     2  stepTime, totalTime, dt, cmname, coordMp, charLength,

     3  props, density, strainInc, relSpinInc,

     4  tempOld, stretchOld, defgradOld, fieldOld,

     5  stressOld, stateOld, enerInternOld, enerInelasOld,

     6  tempNew, stretchNew, defgradNew, fieldNew,

C Write only (modifiable) variables -

     7  stressNew, stateNew, enerInternNew, enerInelasNew )

C

      include 'vaba_param.inc'

C

      dimension props(nprops), density(nblock), coordMp(nblock,*),

     1  charLength(nblock), strainInc(nblock,ndir+nshr),

     2  relSpinInc(nblock,nshr), tempOld(nblock),

     3  stretchOld(nblock,ndir+nshr),

     4  defgradOld(nblock,ndir+nshr+nshr),

     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),

     6  stateOld(nblock,nstatev), enerInternOld(nblock),

     7  enerInelasOld(nblock), tempNew(nblock),

     8  stretchNew(nblock,ndir+nshr),

     8  defgradNew(nblock,ndir+nshr+nshr),

     9  fieldNew(nblock,nfieldv),

     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),

     2  enerInternNew(nblock), enerInelasNew(nblock)

C

      character*80 cmname

C

C  six parameters

      Parameter(EX=2.71828)

      snmax=props(1)  ! the maxinum normal stress

      deltn=props(2)  ! the maxinum normal displacement

      stmax=props(3)   ! the maxinum shear stress

      delts=props(4)  ! the maxinum shear displacement

      t=props(5)   !depth of cohesive elements

      r=props(6)   

 

C

C)

C the maxinum energy

      egyn=EX*snmax*deltn

      egyt=SQRT(EX/2)*stmax*delts

      q=egyt/egyn

C

C)

C

      do 100 i=1, nblock

C

C)))) calculate strain

        sta1=stateOld(i,1)+strainInc(i,1)

        sta2=stateOld(i,2)+strainInc(i,2)

        sta3=stateOld(i,3)+strainInc(i,3)

        sta4=stateOld(i,4)+strainInc(i,4)

C calculate  displacement

        detn=sta1*t

        dets=sta2*t

        

        tn=detn/deltn

        ts=dets/delts

 

C)calculate stress and energy

 

        sig1=(egyn/deltn)*(EX**(-tn))*(tn*(EX**(-ts*ts))+((1-q)/(r-1))

     1   *(1-(EX**(-ts*ts)))*(r-tn))

        sig2=(egyn/deltn)*2*(deltn/delts)*ts*(q+((r-q)/(r-1))*tn)

     2   *(EX**(-tn))*(EX**(-ts*ts))

        sig3=0

        sig4=0

        

        ENGY=egyn+egyn*(EX**(-tn))*((1-r+tn)*((1-q)/(r-1))

     1   -(q+((r-q)/(r-1))*tn)*(EX**(-ts*ts)))

 

C calculate the pure normal and shear energy

        ENt=egyn+egyn*((1-r)*((1-q)/(r-1))

     1   -q*(EX**(-ts*ts)))

        ENn=egyn+egyn*(EX**(-tn))*((1-r+tn)*((1-q)/(r-1))

     1   -(q+((r-q)/(r-1))*tn))

 

C

C................upate the stress ##################################,,

C

        stressNew(i,1)=sig1

        stressNew(i,2)=sig2

        stressNew(i,3)=sig3

        stressNew(i,4)=sig4

        

 

C) calculate the energy difference between the predicted energy and the maxinum energy

        detENt=egyt-ENt

        detENn=egyn-ENn

 

C.....,update the state variables

        stateNew(i,1)=stateOld(i,1)+strainInc(i,1)

        stateNew(i,2)=stateOld(i,2)+strainInc(i,2)

        stateNew(i,3)=stateOld(i,3)+strainInc(i,3)

        stateNew(i,4)=stateOld(i,4)+strainInc(i,4)

        

C   update the state variables

        stateNew(i,5)=2*sta1

        stateNew(i,6)=dets

        stateNew(i,7)=detn

        stateNew(i,8)=ENGY

        stateNew(i,9)=ENt

        stateNew(i,10)=ENn

        stateNew(i,11)=detENt

        stateNew(i,12)=detENn

        

C..........predict the fracture........................,,

        

        if(stateNew(i,11)>0.01) then

          stateNew(i,13)=1

        else

          stateNew(i,13)=0

        endif

 

C...........predict the fracture....................,,

        if(stateNew(i,12)>0.01) then

          stateNew(i,14)=1

        else

          stateNew(i,14)=0

        endif 

 

C..........predict the fracture......,,

        

        if( (stateNew(i,13)==0) .or. (stateNew(i,14)==0) ) then

          stateNew(i,15)=1

        else

          stateNew(i,15)=0

        endif 

        

C

100    continue

C

      return

      end

the abaqus input file is at the attachment