cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c the subroutine "accelerations" computes the accelerations on each particle
c as well as the rate of change of the specific energy (or entropy), u.

        subroutine accelerations (xx,vx,ax,au,hh,p,u,mass,den,divv,
     &     sound,maxmu,n,inorm,alpha,beta,ibulk,alphab,betab,
     &     gamma,ientro,ientro2)

           real * 4 xx(n),vx(n),ax(n),au(n),hh(n),p(n),u(n),mass(n)
           real * 4 den(n),divv(n),sound(n),maxmu(n)

           real * 4 xt(n),vxt(n),axt(n),aut(n),ht(n),pt(n),ut(n)
           real * 4 masst(n),dent(n),divvt(n),soundt(n),maxmut(n)

           real * 4 r(n),v(n),h(n),gradw(n),qi(n),qj(n),cij(n),denij(n)
           real * 4 artv(n),eta(n),mu(n),pforce(n)

           real * 4 alpha,beta,alphab,betab

           ax=0.0      ! initialize the accelerations and du/dt to zero
           au=0.0

           ncount=0    ! a counter

c     store copies of parameters for circular shift

           xt=xx
           vxt=vx
           ht=hh
           axt=ax
           aut=au
           pt=p
           ut=u
           masst=mass
           dent=den
           divvt=divv
           soundt=sound
           maxmu=0.0
           maxmut=0.0

           r=0.0

c     compute forces while any two particles are overlapping

           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))

c     shift the data one location

              xt=cshift(xt,shift=1)
              vxt=cshift(vxt,shift=1)
              axt=cshift(axt,shift=1)
              ht=cshift(ht,shift=1)
              pt=cshift(pt,shift=1)
              ut=cshift(ut,shift=1)
              aut=cshift(aut,shift=1)
              masst=cshift(masst,shift=1)
              dent=cshift(dent,shift=1)
              divvt=cshift(divvt,shift=1)
              soundt=cshift(soundt,shift=1)
              maxmut=cshift(maxmut,shift=1)

              r=xx-xt        ! the distance between the particles

              h=(ht+hh)/2.0  ! their average smoothing length

              v=vx-vxt       ! their relative velocity

c     molecular viscosity

              if (inorm.eq.1) then

                 eta=0.1*(h**2)
                 mu=h*v*r/((r**2)+eta)
                 where ((v*r).gt.0.0) mu=0.0
                 where (abs(mu).gt.maxmu) maxmu=abs(mu)
                 where (abs(mu).gt.maxmut) maxmut=abs(mu)
                 cij=(sound+soundt)/2.0
                 denij=(den+dent)/2.0
                 artv=(-alpha*mu*cij)+(beta*(mu**2))
                 artv=artv/denij
                 
              end if

c     bulk viscosity

              if (ibulk.eq.1) then

                 where (divv.ge.0.0)
                    qi=0.0
                 elsewhere
                     qi=(alphab*hh*den*sound*abs(divv))+
     &                (betab*(hh**2)*den*(divv**2))
                 end where
		 	
                 where (divvt.ge.0.0)
                    qj=0.0
                 elsewhere
                    qj=(alphab*ht*dent*soundt*abs(divvt))+
     &                (betab*(ht**2)*dent*(divvt**2))
                 end where

                 artv=(qi/(den**2))+(qj/(dent**2))

                 where ((r*v).gt.0.0) artv=0.0

              end if

c     next compute the gradient of w, gradw

              call gradww (r,hh,ht,n,gradw)

c     the forces


              if (ientro2.eq.0) then

c              pforce=(p/(den**2))+(pt/(dent**2))    ! the two symmetric forms 
                 pforce=(2.0*sqrt(p*pt)/(den*dent))  ! for the pressure force

                 pforce=(pforce+artv)*gradw          ! add the art. visc. term
                                                     ! and multiply by the 
                                                     ! gradient of the int.
                                                     ! kernal

                 ax=ax-(masst*pforce)                ! add to the running 
                 axt=axt+(mass*pforce)               ! totals
              
              
              else

c     the following is for the case when the entropic formulation of the force
c     is used

                 pforce=ut+((gamma-1.0)*u)
                 pforce=pforce*(den**(gamma-2.0))

                 ax=ax-(masst*(pforce+artv)*gradw)

                 pforce=u+((gamma-1.0)*ut)
                 pforce=pforce*(dent**(gamma-2.0))
                 
                 axt=axt+(mass*(pforce+artv)*gradw)

              end if

c     du/dt

              if (ientro.eq.0) then

                 au=au+(masst*pforce*v)
                 aut=aut+(mass*pforce*v)
              
              else

c     the rate of change of the entropic function

                 au=au+(masst*artv*v*gradw)
                 aut=aut+(mass*artv*v*gradw)
                 
              end if

              ncount=ncount+1
              
           end do

c     shift the temporary totals back to their starting points and add to get
c     the complete sum

           axt=cshift(axt,shift=-ncount)
           aut=cshift(aut,shift=-ncount)

           ax=ax+axt
           au=0.5*(au+aut)

           maxmut=cshift(maxmut,shift=-ncount)
           where (maxmut.gt.maxmu) maxmu=maxmut

        end subroutine accelerations

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "therm_diffusion" computes an additional thermal diffusion
c term which can be added to the energy equation

        subroutine therm_diffusion (xx,vx,hh,mass,den,p,u,divv,sound,n
     &       ,g1,g2,maxmu,gamma,thermdiff)

           real * 4 xx(n),vx(n),hh(n),mass(n),den(n),p(n),u(n)
           real * 4 divv(n),sound(n),thermdiff(n)

           real * 4 xt(n),vxt(n),ht(n),pt(n),ut(n),masst(n),dent(n)
           real * 4 divvt(n),soundt(n),thermdifft(n)
           real * 4 r(n),uu(n),gradw(n),qqi(n),qqj(n)
           real * 4 denij(n),h(n),eta(n),g1,g2

c thermal diffusion

           thermdiff=0.0
           thermdifft=thermdiff
           
           ncount=0
           
           xt=xx
           vxt=vx
           ht=hh
           pt=p
           ut=u
           masst=mass
           dent=den
           divvt=divv
           sound=sqrt(gamma*p/den)
           soundt=sound
           maxmu=0.0

           r=0.0
           
           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))
              
              xt=cshift(xt,shift=1)
              vxt=cshift(vxt,shift=1)
              ht=cshift(ht,shift=1)
              pt=cshift(pt,shift=1)
              ut=cshift(ut,shift=1)
              masst=cshift(masst,shift=1)
              dent=cshift(dent,shift=1)
              divvt=cshift(divvt,shift=1)
              soundt=cshift(soundt,shift=1)
              thermdifft=cshift(thermdifft,shift=1)

              r=xx-xt
              uu=u-ut
              
              qqi=(g1*hh*den*sound)+
     &             (g2*den*(hh**2)*(abs(divv)-divv))
              
              
              qqj=(g1*ht*dent*soundt)+
     &             (g2*dent*(ht**2)*(abs(divvt)-divvt))
              
              qqi=qqi/den
              qqj=qqj/dent
              qqj=(qqi+qqj)/2.0

              denij=(den+dent)/2.0


c     next compute the gradient of w, gradw

              call gradww (r,hh,ht,n,gradw)

              h=(hh+ht)/2.0
              eta=0.1*(h**2)

              where (r.ne.0.0)

              thermdiff=thermdiff+(masst*qqj*uu*r*gradw/
     &             (denij*((r**2)+eta)))
              thermdifft=thermdifft-(mass*qqj*uu*r*gradw/
     &             (denij*((r**2)+eta)))

           end where
                        
           ncount=ncount+1
           
        end do

        thermdifft=cshift(thermdifft,shift=-ncount)
        
        thermdiff=2.0*(thermdiff+thermdifft)
        
        end subroutine therm_diffusion
        
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "push1" advances the particles' positions, velocities, and 
c specific energies forward in time by half a time step

        subroutine push1 (xx,u,vx,ax,au,c1,c2,n)
        
           real * 4 xx(n),u(n),vx(n),ax(n),au(n),c1,c2

           vx=vx+(c1*ax)
           xx=xx+(c2*ax)+(c1*vx)
           u=u+(c1*au)

        end subroutine push1

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "benz_h1" advances the particles' smoothing length forward
c in time by half a time step according to the algorithm of Benz

        subroutine benz_h1 (hh,divv,c1,n)

           real * 4 hh(n),divv(n),c1

           hh=hh+(c1*hh*divv)

        end subroutine benz_h1

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "push2" is called after the accelerations have been computed
c and advances the particles through a complete time step

        subroutine push2 (xx,u,vx,ax,au,c1,c2,delt,xo,vxo,uo,t,n)

           real * 4 xx(n),u(n),vx(n),ax(n),au(n),c1,c2
           real * 4 xo(n),vxo(n),uo(n),delt

           if (t.eq.0.0) then

              xx=xx+(c1*vx)
              vx=vx+(c1*ax)
              u=u+(c1*au)
                        
           else

              vx=vxo+(delt*ax)
              xx=xo+(delt*(vxo+(c1*ax)))
              u=uo+(delt*au)
              
           end if

        end subroutine push2

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "benz_h2" advances the smoothing lengths through a complete
c time step and is called in consort with "push2"

        subroutine benz_h2 (hh,delt,divv,ho,n)

           real * 4 hh(n),divv(n),ho(n)
           real * 4 delt

           if (t.eq.0.0) then

              hh=hh+(c1*divv*hh)

           else

              hh=ho+(delt*divv*hh)

           end if

        end subroutine benz_h2

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "update_den" computes the density according the postions and
c the smoothing lengths of the particles
  
        subroutine update_den (xx,mass,den,hh,n)

           real * 4 xx(n),mass(n),den(n),hh(n)
           real * 4 xt(n),dent(n),masst(n),ht(n)
           
           real * 4 w(n),r(n)

           ncount=0

           xt=xx
           ht=hh
           masst=mass
           den=0.0
           dent=den
        
           r=0.0

           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))

              r=xx-xt   ! distance between particles

              call ww(r,hh,ht,n,w) ! w
              
              den=den+(masst*w) ! running total of the density
              
              if (ncount.gt.0) dent=dent+(mass*w) ! running total of the
                                                  ! shifted density
              
              ncount=ncount+1
              
c     shift the data

              xt=cshift(xt,shift=1)
              ht=cshift(ht,shift=1)
              masst=cshift(masst,shift=1)
              dent=cshift(dent,shift=1)

           end do

           dent=cshift(dent,shift=-ncount)
           den=den+dent
           
        end subroutine update_den

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine update_p computes the pressure for each particle assuming an
c ideal gas

        subroutine update_p (hh,mass,u,den,p,gamma,n,ientro)

           real * 4 hh(n),mass(n),p(n),den(n),u(n),gamma

           if (ientro.eq.0) then
              p=(gamma-1.0)*u*den
           else
              p=u*(den**gamma)
           end if

        end subroutine update_p

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine new_delt computes a new time step

        subroutine new_delt (hh,divv,sound,inorm,alpha,beta,
     &       ibulk,alphab,betab,maxmu,n,c1,c2,cour,delt)

           real * 4 hh(n),divv(n),sound(n),maxmu(n),ht(n),delti(n)
           real * 4 delt,c1,c2,cour
      
c     in the case the molecular viscosity if used

           if (inorm.eq.1.or.(inorm.eq.0.and.ibulk.eq.0)) then

              delti=(hh*abs(divv))+sound+(1.2*((alpha*sound)+
     &             (beta*maxmu)))
              delti=hh/delti
           
           end if

c     in the case the bulk viscosity is used

           if (ibulk.eq.1) then
              
              ht=betab*hh*abs(divv)
              where (divv.ge.0.0) ht=0.0
              delti=(hh*abs(divv))+sound+(1.2*((alphab*sound)+
     &             ht))
              delti=hh/delti
              
           end if
      
           delt=minval(delti)
           delt=cour*delt
           c1=delt/2.0
           c2=(delt**2)/4.0

        end subroutine new_delt

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "reflecting" set the boundry conditions.
c This is done by using 50 guard particles on either end of the domain which 
c are the mirro images of particles in the domain

        subroutine reflecting (xx,vx,xo,vxo,den,mass,p,u,uo,hh,ho,n)

            real * 4 xx(n),vx(n),xo(n),vxo(n),den(n),mass(n),p(n),u(n)
            real * 4 uo(n),hh(n),ho(n)


            xx(1:50)=-xx(101:52:-1)+xx(51)
            vx(1:50)=-vx(101:52:-1)
            xo(1:50)=-xo(101:52:-1)+xo(51)
            vxo(1:50)=-vxo(101:52:-1)
            den(1:50)=den(101:52:-1)
            mass(1:50)=mass(101:52:-1)
            p(1:50)=p(101:52:-1)
            u(1:50)=u(101:52:-1)
            uo(1:50)=uo(101:52:-1)
            hh(1:50)=hh(101:52:-1)
            ho(1:50)=ho(101:52:-1)

            i=n-50
            j=n-52
            k=n-101

            xx(i:n)=xx(n-51)+(xx(n-51)-xx(j:k:-1))
            vx(i:n)=-vx(j:k:-1)
            xo(i:n)=xo(n-51)+(xo(n-51)-xo(j:k:-1))
            vxo(i:n)=-vxo(j:k:-1)
            den(i:n)=-(-den(j:k:-1))
            mass(i:n)=-(-mass(j:k:-1))
            p(i:n)=-(-p(j:k:-1))
            u(i:n)=-(-u(j:k:-1))
            hh(i:n)=-(-hh(j:k:-1))

         end subroutine reflecting

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "double_shock_tube" set up the initial conditions for the
c the double shock tube problem of Woodward and Collela

         subroutine double_shock_tube (xx,vx,p,u,den,mass,xmax,xmin,hh,
     &        partno,gamma,hfac,n,ientro)
         
              real * 4 xx(n),vx(n),p(n),u(n),den(n),mass(n)
              real * 4 hh(n),hfac,xmax,xmin,h

              integer * 4 partno(n)

              h=hfac*(xmax-xmin)/(n-100)

              do i=1,n

                 partno(i)=i
                 hh(i)=h
                 xx(i)=(i-51)*(xmax-xmin)/(n-100)
                 den(i)=1.0
                 mass(i)=den(i)*(xmax-xmin)/(n-100)
                 vx(i)=0.0
                 
                 if (xx(i).lt.0.1) p(i)=1000.0
                 
                 if (xx(i).ge.0.1.and.xx(i).lt.0.9) p(i)=0.01
                 
                 if (xx(i).ge.0.9) p(i)=100.0
                 
                 u(i)=p(i)/den(i)
                 u(i)=u(i)/(gamma-1.0)

                 if (ientro.eq.1) u(i)=p(i)/(den(i)**gamma)

c     here density is mass per unit length

              end do

        end subroutine double_shock_tube
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
c The subroutine "ww" computes the 1-d value of the interpolating kernal
c given a distance (r) between the particles and their smoothing lengths
 
       subroutine  ww (r,hh,ht,n,w)

           real * 4 r(n),hh(n),ht(n),w(n),w2(n),q(n),q1,q2
           
           q1=1.0
           q2=2.0

           q=(abs(r)/hh)

           where (q.le.q1)
        
              w=(2.0/3.0)-(q**2)+(0.5*(q**3))
              w=w/hh

           end where

           where (q.gt.q1.and.q.le.q2)
              
              w=(2.0-q)**3
              w=w/(6.0*hh)

          end where

          where (q.gt.q2) w=0.0

          q=abs(r)/ht
              
          where (q.le.q1) 

             w2=(2.0/3.0)-(q**2)+(0.5*(q**3))
             w2=w2/ht
                 
          end where

          where (q.gt.q1.and.q.le.q2)

              w2=(2.0-q)**3
              w2=w2/(6.0*ht)

           end where

           where (q.gt.q2) w2=0.0

           w=(w+w2)/2.0

        end subroutine ww

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "gradww" computes the gradient of the interpolating kernal

        subroutine gradww (r,hh,ht,n,gradw)

           real * 4 r(n),hh(n),ht(n),gradw(n),gradw2(n),q(n),q1,q2

           q1=1.0
           q2=2.0

           q=(abs(r)/hh)

           where (q.le.q1)

              gradw=(1.5*(r**2)/(hh**3))
              gradw=gradw-(2.0*abs(r)/(hh**2))
              gradw=gradw/hh

           end where

           where (q.gt.q1.and.q.le.q2)

              gradw=(2.0-q)**2
              gradw=-3.0*gradw/hh
              gradw=gradw/(6.0*hh)

           end where 

           where (r.ne.0.0)

              gradw=gradw*r/abs(r)

           elsewhere

              gradw=0.0

           end where

           where (q.gt.q2) gradw=0.0

           q=abs(r)/ht

           where (q.le.q1)

              gradw2=(1.5*(r**2)/(ht**3))
              gradw2=gradw2-(2.0*abs(r)/(ht**2))
              gradw2=gradw2/ht

           end where

           where (q.gt.q1.and.q.le.q2)

              gradw2=(2.0-q)**2
              gradw2=-3.0*gradw2/ht
              gradw2=gradw2/(6.0*ht)

           end where

           where (r.ne.0.0)

              gradw2=gradw2*r/abs(r)

           elsewhere

              gradw2=0.0

           end where

           where (q.gt.q2) gradw2=0.0

           gradw=(gradw+gradw2)/2.0

        end subroutine gradww

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "smooth" smooths a variable (here represented by u).

        subroutine smooth (xx,hh,mass,den,n,u)

           real * 4 xx(n),hh(n),mass(n),den(n),u(n)
           real * 4 xt(n),ht(n),masst(n),old_u(n),old_ut(n),ut(n)
           real * 4 dent(n),r(n),w(n)

           ncount=0

           xt=xx
           ht=hh
           masst=mass
           old_u=u
           old_ut=u
           u=0.0
           ut=0.0

           r=0.0
           
           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))

              r=xx-xt
           
              call ww(r,hh,ht,n,w)
              
              u=u+(masst*old_ut*w)
              
              if (ncount.gt.0) ut=ut+(mass*old_u*w)
              
              ncount=ncount+1
           
              xt=cshift(xt,shift=1)
              ht=cshift(ht,shift=1)
              masst=cshift(masst,shift=1)
              old_ut=cshift(old_ut,shift=1)
              ut=cshift(ut,shift=1)
           
           end do

           ut=cshift(ut,shift=-ncount)
           u=u+ut

           u=u/den

        end subroutine smooth

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "div_v" computes the divergence of the velocity field

        subroutine div_v (xx,vx,hh,mass,den,n,divv)

           real * 4 xx(n),vx(n),hh(n),mass(n),den(n),divv(n)
           real * 4 xt(n),vxt(n),ht(n),divvt(n),masst(n)
           real * 4 r(n),gradw(n),v(n)

           divv=0.0
        
           ncount=0

           xt=xx
           vxt=vx
           ht=hh
           divvt=0.0
           masst=mass

           r=0.0

           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))

              xt=cshift(xt,shift=1)
              vxt=cshift(vxt,shift=1)
              ht=cshift(ht,shift=1)
              divvt=cshift(divvt,shift=1)
              masst=cshift(masst,shift=1)

              r=xx-xt
              v=vx-vxt

c next compute the gradient of w, gradw

              call gradww (r,hh,ht,n,gradw)

              divv=divv+(masst*v*gradw)
c v is in the opposite sense as well as gradw, hence the +
              divvt=divvt+(mass*v*gradw)

              ncount=ncount+1

           end do

           divvt=cshift(divvt,shift=-ncount)
           
           divv=-(divv+divvt)/den

        end subroutine div_v

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine SM_h computes new smoothing lengths according to the alg. of
c Steinmetx and Muller

        subroutine SM_h (xx,hh,mass,den,divv,t,delt,n)

           real * 4 hh(n),xx(n),mass(n),den(n),densmooth(n)
           real * 4 dsmoothdt(n),avh(n),divv(n),aveden

           call smooth_den (xx,hh,mass,den,n,densmooth) ! compute a smoothed 
                                                        ! value of the density
 
           call d_smooth_den (xx,hh,mass,den,densmooth,divv,n,dsmoothdt)

                                                        ! compute the rate of
                                                        ! of change of the 
                                                        ! smoothed density

c 1, estimate the smoothed density

           if (t.eq.0.0) densmooth=den
           densmooth=densmooth+(delt*dsmoothdt)

c 2, calculate the average density of the system

           densmooth=1.0/densmooth
           avden=sum(densmooth)
           avden=1.0/avden
           
c 3, find new smoothing length

           hfac=1.5
           avh=hfac*mass/avden
           
c for 3-d need to use (mass/den)**1/3 and (mass/den)**1/2 for 2-d

           densmooth=1.0/densmooth
           kkk=1.0
           hh=avh*((avden/densmooth)**kkk)

        end subroutine SM_h

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "smooth_den" computes a smoothed value for the density
c for the Steinmetz and Muller alg.

        subroutine smooth_den (xx,hh,mass,den,n,densmooth)

           real * 4 xx(n),hh(n),mass(n),den(n),densmooth(n)
           real * 4 xt(n),ht(n),masst(n),dent(n),denst(n)
           real * 4 r(n),w(n)

c compute a new smoothed density

           ncount=0

           hfac=1.5
           
           xt=xx
           ht=hh
           masst=mass
           dent=den
           densmooth=0.0
           denst=densmooth

           r=0.0

           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))
              
              r=xx-xt

              call ww(r,hh,ht,n,w)
              
              densmooth=densmooth+(masst*w*dent)
              
              if (ncount.gt.0) denst=denst+(mass*w*den)
              
              ncount=ncount+1
              
              xt=cshift(xt,shift=1)
              ht=cshift(ht,shift=1)
              masst=cshift(masst,shift=1)
              dent=cshift(dent,shift=1)
              denst=cshift(denst,shift=1)
              
           end do

           denst=cshift(denst,shift=-ncount)
           densmooth=(densmooth+denst)/den

           hh=hfac*(mass/densmooth)

        end subroutine smooth_den

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c The subroutine "d_smooth_den" computes the rate of change of the smoothed
c density for the Steinmetz + Muller alg.

        subroutine d_smooth_den (xx,hh,mass,den,densmooth,divv,n,
     &       dsmoothdt)

           real * 4 xx(n),hh(n),mass(n),den(n),divv(n),dsmoothdt(n)
           real * 4 xt(n),ht(n),masst(n),dent(n),divvt(n),dsmt(n)
           real * 4 w(n),r(n),densmooth(n)

           ncount=0

           xt=xx
           ht=hh
           masst=mass
           dent=densmooth
           divvt=divv
           dsmoothdt=0.0
           dsmt=0.0

           r=0.0
           
           do while (any(abs(r).le.(2.0*hh).or.abs(r).le.(2.0*ht)))
           
              r=xx-xt
              
              call ww(r,hh,ht,n,w)
              
              dsmoothdt=dsmoothdt+(masst*dent*divvt*w)
              
              if (ncount.gt.0) dsmt=dsmt+(mass*densmooth*divv*w)

              ncount=ncount+1
              
              xt=cshift(xt,shift=1)
              ht=cshift(ht,shift=1)
              masst=cshift(masst,shift=1)
              dent=cshift(dent,shift=1)
              divvt=cshift(divvt,shift=1)
              dsmt=cshift(dsmt,shift=1)

           end do

           dsmt=cshift(dsmt,shift=-ncount)
           dsmoothdt=(dsmoothdt+dsmt)
           
           dsmoothdt=2.0*dsmoothdt/den
           dsmoothdt=(densmooth*divv)-dsmoothdt

        end subroutine d_smooth_den
           
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc