changing density for FGM in umat

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