40       subroutine bardmc (y,      shgl,      shp,   argument
71      &          shgl(MAXTOP,nsd,maxsh,MAXQPT), shp(MAXTOP,maxsh,MAXQPT),    
251      &               shp(lcsyst,1:nshl,:),Qwtf(lcsyst,1:ngaussf))
270      &               shp(lcsyst,1:nshl,:))
561      &        mxmudmi(iblk)%p,cdelsq,shp(lcsyst,1:nshl,:))
672       subroutine hfilter (y, x, ien, hfres, shgl, shp, Qwtf)  argument
684      &          shp(nshl,maxsh),           Qwtf(ngaussf),
756             rho = rho+shp(i,intp)*yl(:,i,1) !density at qpt
765             u1 = u1 + shp(i,intp)*yl(:,i,2)
766             u2 = u2 + shp(i,intp)*yl(:,i,3)
767             u3 = u3 + shp(i,intp)*yl(:,i,4)
773             fresli(:,i,1) = rho * u1 * shp(i,intp) !rho * u1 * WdetJ
774             fresli(:,i,2) = rho * u2 * shp(i,intp) !rho * u2 * WdetJ
775             fresli(:,i,3) = rho * u3 * shp(i,intp) !rho * u3 * WdetJ
776             fresli(:,i,4) = rho * shp(i,intp) !rho * WdetJ
777             fresli(:,i,5) = WdetJ*shp(i,intp) !Integral of filter kernel
779             fresli(:,i,6) = rho * u1 * u1 * shp(i,intp) !rho * u1 * u1 * WdetJ
780             fresli(:,i,7) = rho * u2 * u2 * shp(i,intp) !rho * u2 * u2 * WdetJ
781             fresli(:,i,8) = rho * u3 * u3 * shp(i,intp) !rho * u3 * u3 * WdetJ
782             fresli(:,i,9) = rho * u1 * u2 * shp(i,intp) !rho * u1 * u2 * WdetJ
783             fresli(:,i,10)= rho * u1 * u3 * shp(i,intp) !rho * u1 * u3 * WdetJ
784             fresli(:,i,11)= rho * u2 * u3 * shp(i,intp) !rho * u2 * u3 * WdetJ
826      &     hfres, shgl, shp, Qwtf)  argument
838      &          shp(nshl,maxsh),
914         fresli(:,22) = fresli(:,22)+shp(i,intp)*yl(:,i,1)     ! density at qpt
915         fresli(:,25) = fresli(:,25)+shp(i,intp)*hfresl(:,i,1) !bar(rho u1)
916         fresli(:,26) = fresli(:,26)+shp(i,intp)*hfresl(:,i,2) !bar(rho u2)
917         fresli(:,27) = fresli(:,27)+shp(i,intp)*hfresl(:,i,3) !bar(rho u3)
1003        u1 = u1 + shp(i,intp)*yl(:,i,2)
1004        u2 = u2 + shp(i,intp)*yl(:,i,3)
1005        u3 = u3 + shp(i,intp)*yl(:,i,4)
1060       subroutine disPtion (y, x, ien, shgl, shp, eps, xmudmi)   argument
1069      &          shgl(nsd,nshl,maxsh),       shp(nshl,maxsh),
1161         fresli(:,22) = fresli(:,22)+shp(i,intp)    ! unit density at qpt