changing density for FGM in umat

[sdata]

Hello friends,
Here is the Abaqus user subroutine for Functionally Graded Material with direct integration method. This Umat is based on S.H.Song's Umat. Obviously the coefficient AAA in calculating E is the ratio of young moduli of two sides of the plate. With this method, one can change young modulus through the thickness.

But I have problem with changing density through the thickness, meaning from density of ceramic to density of metal.
I would be grateful if anyone can help me. Any help would be appreciated.

*USER SUBROUTINES
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATEV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATEV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
C
PARAMETER (ONE=1.0D0, TWO=2.0D0)
C Determine material properties based on global coordinates of gauss points.
C COORDS(1) is X-coordinates of gauss points.
C COORDS(2) is Y-coordinates of gauss points.
C COORDS(3) is Z-coordinates of gauss points.
C PROPS is defined by users.
C The function can be also defined by users.
C
C------------- CHNGE -------------------------------------------------------------------------------C
C---------------------------------------------------------------------------------------------------C
E=PROPS(1)*EXP(AAA*COORDS(1))
V=PROPS(2)
C---------------------------------------------------------------------------------------------------C
C---------------------------------------------------------------------------------------------------C
C
C Determine Lame's constantes
C amu is mu
C alambeda is lambda
AMU=E/2.0d0/(1.0d0+V)
ALAMBDA=E*V/(1.0d0+V)/(1.0d0-2*V)
C
C
C Determine the tangant or Jacobian matrix
C
C NDI is number of normal stresses
C NSHR is number of shear stresses
C NTENS is dimension of constitutive matrix
C
C

DO I=1,NDI
DO J=1,NDI
if (I.eq.J) then
ddsdde(I,I)=ALAMBDA+2*AMU
else
ddsdde(I,J)=ALAMBDA
endif
ENDDO
ENDDO
DO I=NDI+1,NTENS
ddsdde(I,I)=AMU
ENDDO
C
C Determine the stress and update the stress
C
DO I=1,NTENS
DO J=1,NTENS
STRESS(I)=STRESS(I)+ddsdde(I,J)*DSTRAN(J)
ENDDO
ENDDO
C
RETURN
END