关于marc的一个子程序问题

麻烦大家看一我用这个写的损伤值,但是结果出来好像不对,麻烦大家看一下!!!
subroutine ueldam(m,n,nn,kc,inc,lovl,matus,timinc,cptim,toten,
2 deven,totend,totenv,surfc,surfd,dt,dtdl,
3 damd,ddamd)
c
include '../common/implicit'
dimension dt(*),dtdl(*),matus(2)
c
c
c User subroutine to determine value of damage parameter
c
c Input:
c
c m = element number
c n = element/elsto number
c nn = integration point number
c kc = layer number
c inc = increment number
c lovl = 4 for assembly phase, 6 for stress recovery phase
c matus(1) = user material set number
c matus(2) = internal material set number
c timinc = time increment
c cptim = time at beginning of increment
c toten = total instantaneous strain energy at the end
c of the current step excluding damage
c deven = deviatoric part of the instantaneous strain energy at
c the end of the current step excluding damage
c totend = stored deviatoric energy at previous step (including damage)
c totenv = stored volumetric energy at previous step (including damage)
c surfc = current radius of continuous damage surface
c surfd = current radius of discontinuous damage surface
c dt = temperature
c dtdl = incremental temperature
c
c Required Output:
c
c damd = value of Kachanov deviatoric damage parameter
c ddamd = derivative of damage parameter with respect to maximum
c total strain energy
c
c
c include 'concom'
common/rubber/WW,vvv,fovinc(1600,9),fov(1600,9)
common/rubber/endom(1600,9)
vvv=toten
dinfalf=0.7719389
nualf=1.2
kc=1.D0
write(6,*) 'dans la routine'
call elmvar (1,m,nn,kc,strain)
fov(m,nn)=strain+1
if(abs(fovinc(m,nn)).gt.abs(fov(m,nn))) then
fov(m,nn)=fovinc(m,nn)
else
fov(m,nn)=fov(m,nn)
endif
damd=1-dinfalf*(1-exp((1-fov(m,nn))*nualf))
endom(m,nn)=damd
return
end