#include "cmhog.def"
c=======================================================================
c/////////////////////////  SUBROUTINE HDFALL  \\\\\\\\\\\\\\\\\\\\\\\\\
c
      subroutine hdfall(filename)
c
c  MAKES AN HDF DUMP CONTAINING ALL ACTIVE FIELD ARRAYS
c
c     written by: Jim Stone
c     date:       January,1989
c     modified1:  PJT - april 1995: 3d -> 2d for the isothermal bar project
c
c  PURPOSE: Makes an hdf dump of all the active field variables.
c    Data is written in the Scientific Data Set to "filename".  Note
c    that data must be stored contiguously in order to interface
c    correctly to the C hdf routines.
c
c  EXTERNALS: HDF library routines, PGAS (only if non-isothermal)
c
c  LOCALS:
c-----------------------------------------------------------------------
c      implicit NONE
#include "param.h"
#include "grid.h"
#include "field.h"
#include "root.h"
      character*8  filename
c
      integer  i,j,k,rank,shape(3),ret,nyd,nzd
      real data(jn*kn),yscale(jn),zscale(kn)
      character*32 string
      equivalence (data,wijk0)
c
      integer  dssdims,dssdast,dssdisc,dsadata,dspdata
      external dssdims,dssdast,dssdisc,dsadata,dspdata
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////
c=======================================================================
c
      do 20 j=jbtm,jtp
        yscale(j-jbtm+1) = y(j) + 0.5*dy(j)
20    continue
      do 30 k=kbtm,ktp
        zscale(k-kbtm+1) = z(k) + 0.5*dz(k)
30    continue
      nyd = jtp-jbtm+1
      nzd = ktp-kbtm+1
c
c  double check, since the two versions (stone vs. piner) seem to use 
c  different indices
c
      if (nyd .ne. je-js+1) then
         write(*,*) 'Problem with Radial grid',nyd,js,je         
      endif
      if (nzd .ne. ke-ks+1) then
         write(*,*) 'Problem with Angular grid',nzd,ks,ke
      endif


c      rank     = 3
      rank     = 2
      shape(1) = nyd
      shape(2) = nzd
      ret = dssdims(rank,shape)
c      ret = dssdisc(1,shape(1),xscale)
c      ret = dssdisc(2,shape(2),yscale)
c      ret = dssdisc(3,shape(3),zscale)
      ret = dssdisc(1,shape(1),yscale)
      ret = dssdisc(2,shape(2),zscale)

c
c  y-velocity
c
      do 200 k=ks,ke
      do 200 j=js,je
        data((k-ks)*nyd + (j-js)+1) = v(j,k)
200   continue   
      write(string,"('R-VELOCITY AT TIME=',1pe8.2)") time
      ret = dssdast(string,'km/sec   ','1pe8.2','Polar')
      ret = dsadata(filename,rank,shape,data)
c
c  z-velocity
c
      do 300 k=ks,ke
      do 300 j=js,je
        data((k-ks)*nyd + (j-js)+1) = w(j,k)
300   continue   
      write(string,"('PHI-VELOCITY AT TIME=',1pe8.2)") time
      ret = dssdast(string,'km/sec   ','1pe8.2','Polar')
      ret = dsadata(filename,rank,shape,data)
c
c  density
c
      do 400 k=ks,ke
      do 400 j=js,je
        data((k-ks)*nyd + (j-js)+1) = d(j,k)
400   continue
      write(string,"('DENSITY AT TIME=',1pe8.2)") time
      ret = dssdast(string,'Msolar/pc**2','1pe8.2','Polar')
      ret = dsadata(filename,rank,shape,data)

      return
      end