Several questions for VUMAT

Hi all, I am a new ABAQUS VUMAT user. Recently I try to define the matreial behavior in stamping of sheet metal  through the VUMAT. there are several questions about VUMAT.

（1）accoring to the ABAQUS user manual，for plane stress problem，stressnew(nblock,ndir+nshr), where ndir=3,nshr=1, but in fact,there are 5 stress commponents: sigma_{11},sigma_{22},sigma_{12},sigma_{23},sigma_{31}. how to match stressnew(*,i)(i=1,2,3,4) with the stress commponents: sigma_{11},sigma_{22},sigma_{12},sigma_{23},sigma_{31}?

in my view, the direct commponents：tressnew(*,1)=sigma_{11},stressnew(*,2)=sigma_{22},stressnew(*,3)=sigma_{33} and the indirect commponent：stressnew(*,4)=sigma_{12}. However,how about sigma_{23}and sigma_{31}?。

how to update the  sigma_{23}and sigma_{31}?

Is it right to use the transverse shear stiffness update thesigma_{23}and sigma_{31} ? for example K*（epsilon_{23},epsilon_{31}） where K means the transverse shear stiffness？
（2）the following is a simple example of the coding of subroutine VUMAT in ABAQUS user manual. there are two questions:
where was STATE(*,i)(i=1,2,3,4,5) defined?andwhere was props(i)(i=1,2,3,4) defined？
subroutine vumat(
C Read only –
     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,
     3  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only –只写
     5  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
C J2 Mises Plasticity with kinematic hardening for plane
C strain case.
C Elastic predictor, radial corrector algorithm.
C
C The state variables are stored as:
C      STATE(*,1) = back stress component 11 ***（where was STATE(*,1) defined?）
C      STATE(*,2) = back stress component 22 ****（where was STATE(*,2) defined？）
C      STATE(*,3) = back stress component 33 ****（where was STATE(*,3) defined？）
C      STATE(*,4) = back stress component 12 ****（where was STATE(*,4) defined？）
C      STATE(*,5) = equivalent plastic strain *****（where was STATE(*,5) defined？）
C
C
C All arrays dimensioned by (*) are not used in this algorithm
      dimension props(nprops), density(nblock),
     1  coordMp(nblock,*),
     2  charLength(*), strainInc(nblock,ndir+nshr),
     3  relSpinInc(*), tempOld(*),
     4  stretchOld(*), defgradOld(*),
     5  fieldOld(*), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(*),
     8  stretchNew(*), defgradNew(*), fieldNew(*),
     9  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     1  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
      parameter( zero = 0., one = 1., two = 2., three = 3.,
     1  third = one/three, half = .5, twoThirds = two/three,
     2  threeHalfs = 1.5 )
C
      e     = props(1)*****（where was props(1) defined？）
      xnu   = props(2)*****（where was props(2) defined？）
      yield = props(3)*****（where was props(3) defined？）
      hard  = props(4)*****（where was props(4) defined？）

C
      twomu  = e / ( one + xnu )
      thremu = threeHalfs * twomu
      sixmu  = three * twomu
      alamda = twomu * ( e - twomu ) / ( sixmu - two * e )
      term   = one / ( twomu * ( one + hard/thremu ) )
      con1   = sqrt( twoThirds )
C
      do 100 i = 1,nblock
C
C Trial stress
        trace  = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
        sig1 = stressOld(i,1) + alamda*trace + twomu*strainInc(i,1)
        sig2 = stressOld(i,2) + alamda*trace + twomu*strainInc(i,2)
        sig3 = stressOld(i,3) + alamda*trace + twomu*strainInc(i,3)
        sig4 = stressOld(i,4)                + twomu*strainInc(i,4)
C
C Trial stress measured from the back stress
        s1 = sig1 - stateOld(i,1)
        s2 = sig2 - stateOld(i,2)
        s3 = sig3 - stateOld(i,3)
        s4 = sig4 - stateOld(i,4)
C
C Deviatoric part of trial stress measured from the back stress
        smean = third * ( s1 + s2 + s3 )
        ds1 = s1 - smean
        ds2 = s2 - smean
        ds3 = s3 - smean
C
C Magnitude of the deviatoric trial stress difference
        dsmag = sqrt( ds1**2 + ds2**2 + ds3**2 + 2.*s4**2 )
C
C Check for yield by determining the factor for plasticity,
C zero for elastic, one for yield
        radius = con1 * yield
        facyld = zero
        if(  dsmag - radius .ge. zero ) facyld = one
C
C Add a protective addition factor to prevent a divide by zero
C when dsmag is zero. If dsmag is zero, we will not have exceeded
C the yield stress and facyld will be zero.
        dsmag  = dsmag + ( one - facyld )
C
C Calculated increment in gamma (this explicitly includes the
C time step)
        diff   = dsmag - radius
        dgamma = facyld * term * diff
C
C Update equivalent plastic strain
        deqps  = con1 * dgamma
        stateNew(i,5) = stateOld(i,5) + deqps
C
C Divide dgamma by dsmag so that the deviatoric stresses are
C explicitly converted to tensors of unit magnitude in the
C following calculations
        dgamma = dgamma / dsmag
C
C Update back stress
        factor  = hard * dgamma * twoThirds
        stateNew(i,1) = stateOld(i,1) + factor * ds1
        stateNew(i,2) = stateOld(i,2) + factor * ds2
        stateNew(i,3) = stateOld(i,3) + factor * ds3
        stateNew(i,4) = stateOld(i,4) + factor *  s4
C
C Update the stress
        factor   = twomu * dgamma
        stressNew(i,1) = sig1 - factor * ds1
        stressNew(i,2) = sig2 - factor * ds2
        stressNew(i,3) = sig3 - factor * ds3
        stressNew(i,4) = sig4 - factor *  s4
C
C Update the specific internal energy -
        stressPower = half * (
     1    ( stressOld(i,1)+stressNew(i,1) )*strainInc(i,1)
     1    +     ( stressOld(i,2)+stressNew(i,2) )*strainInc(i,2)
     1    +     ( stressOld(i,3)+stressNew(i,3) )*strainInc(i,3)
     1    + two*( stressOld(i,4)+stressNew(i,4) )*strainInc(i,4) )
C
        enerInternNew(i) = enerInternOld(i)
     1    + stressPower / density(i)
C
C Update the dissipated inelastic specific energy -
        smean = third * ( stressNew(i,1) + stressNew(i,2)
     1                + stressNew(i,3) )
        equivStress = sqrt( threeHalfs *
     1    ( (stressNew(i,1)-smean)**2
     1    + (stressNew(i,2)-smean)**2
     1    + (stressNew(i,3)-smean)**2
     1    + two * stressNew(i,4)**2 ) )
        plasticWorkInc = equivStress * deqps
        enerInelasNew(i) = enerInelasOld(i)
      1   + plasticWorkInc / density(i)
  100 continue
C
      return
      end