!  $Id: hd.f,v 1.12 2000/09/24 19:44:09 aake Exp $
!***********************************************************************
      subroutine pde
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      include 'stagger.ifc'

      character*64:: id=
     & "$Id: hd.f,v 1.12 2000/09/24 19:44:09 aake Exp $"

      if (id .ne. '') then
      print *,id
      id = ''
      end if

      cpu = elapsed()

!-----------------------------------------------------------------------
!  Face centered velocities, boundary conditions
!-----------------------------------------------------------------------
!     call extrap_full_log (r)
!     call periodic(r)
      call mixed_full(r,ssbot)
!     call extrap_full_log (e)
!     call periodic(e)
      call mixed_full(e,eebot)
!     call antisymm_half (px)
!     call symmetric_half (px)
!     call extrap_half (px)
!     call periodic(px)
!     call set_half(Ux,Ux(lb+1,:,:),eebot)
      xdnr = exp(xdn(alog(r)))
      ydnr = ydn(r)
      zdnr = zdn(r)
      Ux = px/xdnr
      Uy = py/ydnr
      Uz = pz/zdnr
!     call symmetric_half (Ux)
      call mixed_half(e,eebot)
  	px = Ux*xdnr
      call symmetric_full (Uy) 
      call symmetric_full (Uz)
      nflop=nflop+2*5+2

!-----------------------------------------------------------------------
!  Pressure (pp), thermal energy per unit mass (ee) and sound speed (cs)
!-----------------------------------------------------------------------
!     e(ub,:,:) = r(ub,:,:)*eebot
      ee = e/r
      pp = e*(gamma-1.)
      nflop=nflop+5+4

!-----------------------------------------------------------------------
!  Numerical diffusion coefficients.  With Sij = ui,j + uj,i good values
!  are [nu1,nu2]=[0.03,0.03] ([0.02,0.02] and [0.01,0.03] are too small) 
!  and nu3=0.3-0.4 (0.2 is too small).
!-----------------------------------------------------------------------
      cs = sqrt((gamma*pp)/r)
      nflop=nflop+12

      scr1 = amax1((difx1(Ux)+dify1(Uy)+difz1(Uz)),0.0)
      scr1 = smooth3(max3(scr1))
      call symmetric_full (scr1)

      vx = (nu1*dx)*xdn1(cs)+(nu2*dx)*abs(Ux)+(nu3*dx)*xdn1(scr1)
      vy = (nu1*dy)*ydn1(cs)+(nu2*dy)*abs(Uy)+(nu3*dy)*ydn1(scr1)
      vz = (nu1*dz)*zdn1(cs)+(nu2*dz)*abs(Uz)+(nu3*dz)*zdn1(scr1)
      nflop=nflop+18

!-----------------------------------------------------------------------
!  Mass conservation, cell centered
!-----------------------------------------------------------------------
      drdt =  drdt
     &        - ddxup(px)
     &        - ddyup(py) 
     &        - ddzup(pz)
      nflop=nflop+3

!-----------------------------------------------------------------------
!  Equations of motion, XX
!-----------------------------------------------------------------------
      S = ddxup(Ux)*2.
      TT = r*xup1(vx)*S - pp
      dEdt = dEdt + TT*S*0.5
      TT = xup(Ux)**2*r - TT
      dpxdt = dpxdt - ddxdn(TT) + grav*xdnr
      nflop=nflop+9

!-----------------------------------------------------------------------
!  optional radial mode damping
!-----------------------------------------------------------------------
      if (pxdamp .gt. 0.) then
          do i=1,mx
            pxav = sum(px(i,:,:))/(my*mz)
            dpxdt(i,:,:) = dpxdt(i,:,:) - pxav/pxdamp
          enddo
          nflop=nflop+5
      endif

!-----------------------------------------------------------------------
!  YY
!-----------------------------------------------------------------------
      S = ddyup(Uy)*2.
      TT = r*yup1(vy)*S - pp
      dEdt = dEdt + TT*S*0.5
      TT = yup(Uy)**2*r - TT
      dpydt = dpydt - ddydn(TT) 
      nflop=nflop+9
 
!-----------------------------------------------------------------------
!  ZZ
!-----------------------------------------------------------------------
      S = ddzup(Uz)*2.
      TT = r*zup1(vz)*S - pp
      dEdt = dEdt + TT*S*0.5
      TT = zup(Uz)**2*r  - TT
      dpzdt = dpzdt - ddzdn(TT)
      nflop=nflop+9
 
!-----------------------------------------------------------------------
!  XY
!-----------------------------------------------------------------------
      S = ddydn(Ux)+ddxdn(Uy)
      TT = ydn(xdnr)*0.5*(ydn1(vx)+xdn1(vy))*S
      dEdt = dEdt + xup1(yup1(TT*S))
      TT = ydn(px)*xdn(Uy) - TT
      dpxdt = dpxdt - ddyup (TT)
      dpydt = dpydt - ddxup (TT)
      nflop=nflop+13

!-----------------------------------------------------------------------
!  XZ
!-----------------------------------------------------------------------
      S = ddxdn(Uz)+ddzdn(Ux)
      TT = zdn(xdnr)*0.5*(zdn1(vx)+xdn1(vz))*S
      dEdt = dEdt + xup1(zup1(TT*S))
      TT = zdn(px)*xdn(Uz) - TT
      dpxdt = dpxdt - ddzup (TT)
      dpzdt = dpzdt - ddxup (TT)
      nflop=nflop+13

!-----------------------------------------------------------------------
!  YZ
!-----------------------------------------------------------------------
      S = ddydn(Uz)+ddzdn(Uy)
      TT = ydn(zdnr)*0.5*(zdn1(vy)+ydn1(vz))*S
      dEdt = dEdt + yup1(zup1(TT*S))
      TT = zdn(py)*ydn(Uz) - TT
      dpydt = dpydt - ddzup (TT)
      dpzdt = dpzdt - ddyup (TT)
      nflop=nflop+13

!-----------------------------------------------------------------------
!  Energy flux divergence, cell centered
!-----------------------------------------------------------------------
      if (tcool .eq. 0.) then
          Qcool = 0.
!         e(lb,:,:) = eetop*r(lb,:,:)
!         e(ub,:,:) = eebot*r(ub,:,:)
      else
        do i=1,5
            Qcool(i,:,:) = r(i,:,:)*(ee(i,:,:)-eetop)*fcool(i)*(7-i)
        enddo
        do i=6,mx/2-1
            Qcool(i,:,:) = r(i,:,:)*(ee(i,:,:)-eetop)*fcool(i)
        enddo
        do i=mx/2,mx
            Qcool(i,:,:) = r(i,:,:)*(ee(i,:,:)-x(i))*fcool(i)
        enddo
          nflop=nflop+3
      endif
      dEdt = dEdt + Qcool
     &      - ddxup (px*xdn(ee) - xdnr*vx*ddxdn(ee))
     &      - ddyup (py*ydn(ee) - ydnr*vy*ddydn(ee)) 
     &      - ddzup (pz*zdn(ee) - zdnr*vz*ddzdn(ee))
      nflop=nflop+16

      zones = zones + mx*my*mz
      cpu = elapsed()
      sec = sec + cpu

      END
