subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,props        
     1 		,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,kstep,
     2 		kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,lflags
     3 		,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)
c
		include 'aba_param.inc'
c
		dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),svars(nsvars),
     1 		energy(8),props(*),coords(mcrd,nnode),
     2 		u(ndofel),du(mlvarx,*),v(ndofel),a(ndofel),time(2),
     3 		params(3),jdltyp(mdload,*),adlmag(mdload,*),
     4 		ddlmag(mdload,*),predef(2,npredf,nnode),lflags(*),jprops(*)
c
		dimension gx(4),hx(4),phi(4),phix(4),phiy(4),phig(4),phih(4)
		dimension rjac(2,2),rjaci(2,2)
c
c     	material property definition
		rho   = 1.
		spec  = 1.
		cond = 1.
c     	initialization  (nrhs=1)                                                         
		do k1=1,ndofel                      
			rhs(k1,nrhs)=0.
			do k2=1,ndofel
				amatrx(k2,k1)=0.
			enddo                                    
		enddo                                      
		if (lflags(3).eq.4) return
c     	transient analysis
		if (lflags(1).eq.33) then
c     		determine gauss point locations
			gpos=1./sqrt(3.)	
			gx(1)=-gpos
			gx(2)=gpos
			gx(3)=gpos
			gx(4)=-gpos
			hx(1)=-gpos
			hx(2)=-gpos
			hx(3)=gpos
			hx(4)=gpos
c     		assemble amatrx and rhs
			do k=1,4
c    			loop through gauss pts
				g=gx(k)
				h=hx(k)
c				get shape functions and derivatives			
				phi(1)=0.25*(1.-g)*(1.-h)
				phi(2)=0.25*(1.+g)*(1.-h)
				phi(3)=0.25*(1.+g)*(1.+h)
				phi(4)=0.25*(1.-g)*(1.+h)
				phig(1)=0.25*-(1.-h)
				phig(2)=0.25*(1.-h)
				phig(3)=0.25*(1.+h)
				phig(4)=0.25*-(1.+h)
				phih(1)=0.25*-(1.-g)
				phih(2)=0.25*-(1.+g)
				phih(3)=0.25*(1.+g)
				phih(4)=0.25*(1.-g)				
c				get ip coords
				crdx=0.
				crdy=0.
				do k1=1,nnode
					crdx=crdx+phi(k1)*coords(1,k1)
					crdy=crdy+phi(k1)*coords(2,k1)
				end do  
c				get jacobian				
				rjac=0.
				do i=1,4
					rjac(1,1)=rjac(1,1)+phig(i)*coords(1,i)
					rjac(1,2)=rjac(1,2)+phig(i)*coords(2,i)
					rjac(2,1)=rjac(2,1)+phih(i)*coords(1,i)
					rjac(2,2)=rjac(2,2)+phih(i)*coords(2,i)
				enddo
				djac=rjac(1,1)*rjac(2,2)-rjac(1,2)*rjac(2,1)
				print *,djac
				rjaci(1,1) =  rjac(2,2)/djac
				rjaci(2,2) =  rjac(1,1)/djac
				rjaci(1,2) = -rjac(1,2)/djac
				rjaci(2,1) = -rjac(2,1)/djac
c				get b matrix
				phix(1)=rjaci(1,1)*phig(1)+rjaci(1,2)*phih(1)
				phiy(1)=rjaci(2,1)*phig(1)+rjaci(2,2)*phih(1)
				phix(2)=rjaci(1,1)*phig(2)+rjaci(1,2)*phih(2)
				phiy(2)=rjaci(2,1)*phig(2)+rjaci(2,2)*phih(2)
				phix(3)=rjaci(1,1)*phig(3)+rjaci(1,2)*phih(3)
				phiy(3)=rjaci(2,1)*phig(3)+rjaci(2,2)*phih(3)
				phix(4)=rjaci(1,1)*phig(4)+rjaci(1,2)*phih(4)
				phiy(4)=rjaci(2,1)*phig(4)+rjaci(2,2)*phih(4)																	
c				interpolate temperatures to int points
				dtdx=u(1)*phix(1)+u(2)*phix(2)
     1				+u(3)*phix(3)+u(4)*phix(4)
				dtdy=u(1)*phiy(1)+u(2)*phiy(2)
     1				+u(3)*phiy(3)+u(4)*phiy(4)	
				t=u(1)*phi(1)+u(2)*phi(2)+u(3)*phi(3)+u(4)*phi(4)
				told=(u(1)-du(1,nrhs))*phi(1)+(u(2)-du(2,nrhs))*phi(2)+
     1				(u(3)-du(3,nrhs))*phi(3)+(u(4)-du(4,nrhs))*phi(4)
c				other housekeeping			
				dtdt=(t-told)/dtime				
				we=djac
c				Assemble Element Stiffness Matrix and Add to Global
				do ki=1,4
c       			loop over nodes
					rhs(ki,nrhs)=rhs(ki,nrhs)-we*(phi(ki)*rho*spec*dtdt 
     1  				+ cond*(phix(ki)*dtdx+phiy(ki)*dtdy))
					do kj=1,4
						amatrx(ki,kj)=amatrx(ki,kj)+we*(phi(ki)*phi(kj)
     1          			*rho*spec/dtime+cond*(phix(ki)*phix(kj)+
     1 						phiy(ki)*phiy(kj)))
					end do
				end do
			enddo
		end if
		return
		end