! $Id: unf_to_hdf.f90,v 1.24 2000/11/21 22:14:06 crumleyj Exp $

      module unf_to_hdf_mod

      contains

C This subroutine pulls data from an unformatted file and writes it to
C an hdf file.  It splits the unformatted dat into sections for each iteration.
C The unformatted binary file will (at least at this point) only carry the
C actual numerical data.  Other useful data items such as labels for 
c variables, and titles for data sets will be included in the HDF file.
c This information will be reconstructed in this routine, based on the 
c what's input into this subroutine.

      subroutine unf_to_hdf (filename, name_len, unf_unit, type,
     #    data_name, dat_name_len, num_sets,it_0,it_stride, average,
     #    ave_over, drop_first_set, sd_id_input) 

C    filename is the name of the hdf file, unf_unit is the unit of the
c    unformatted file, type is the data type of the file and size of the
C    data array. Data_name is the name of the data, it will be combined
c    with other info to make the set label.  Num_sets, it_0, and it_stride
c    are the number of sets of data in a file, it_0 is the first iteration
c    that a data set is saved for, it_stride is the number of iterations
c    between 2 consecutive data sets.  Average is an optional logical parameter
c    that is set by the calling subroutine if unf_to_hdf is supposed to
c    average the data down to make a smaller hdf file.  This parameter
c    is OPTIONAL. Ave_over is also optional.  It is the number of iterations
c    to average over.  Drop_first_set is an optional parameter that tells 
c    the subroutine to read in and drop the first data set.  This is for 
C    imitating the behavior of output - the first phi data set (it 0) is not
C    in making phiav.  In_filename is an optional parameter that specifies the
c    input hdf filename to read from - only used when averaging data down.

C This subroutine pulls alot from wr_xory_hdf.f and the other wr*hdf*.f's.
C  
      
C     Include file that contains hdf.h and other definitions
      include 'maxpar.h' 
      include 'param.h'
      include 'hist.h'
      include 'hdf.inc'
c      include  '/net/belka/d3/crumley/local/hdf/hdf/include/hdf.inc'

      integer*4 :: sd_id, sds_id, rank, unf_unit, type, ios_read_unf=0
      integer, parameter :: REAL_HDF_KIND=DFNT_FLOAT32, 
     # 		REAL_OUTPUT_KIND = 4, DEFAULT_KIND=8
c  REAL_HDF_KIND is a parameter that will be used with HDF dataset calls.
c  The values that it can be set at are DFNT_FLOAT32 or DFNT_FLOAT64.  It
C  tells hdf how big the data data quantities are.
C  REAL_OUTPUT_KIND is the number of bytes in the real output.  It is used
C  to set the size of output arrays used with hdf.

      integer sfstart, sfcreate, sfendacc, sfend, dimsizes(2)
      integer status,dim_id1,dim_id0
      integer sfsdmstr,sfdimid
      parameter (ngmax=ngymax)
      character*(maxchar) set_label
      
      integer :: itnum=0,i=0,jj=0
      character *(maxchar) filename,setlabel, data_name
      integer start(2),edges(2),stride(2),sfsdmname
      integer dat_name_len,it_0, it_stride,it_num,it_num_len
      character*(maxdigits) :: it_num_ch, spec_char

      !Below is an array that holds the data sizes. The funny reshape
      !command is used instead of a data statement, to initialize the
      !array.  Values are listed so that the ith and ith + num_type_file
      !values are the sizes for file i. (Use num_files -1, because 3-d
      ! data not considered here - so its not really all sizes anymore).
      integer,dimension(num_type_file-1,2) :: all_sizes=reshape
     #    ((/2,ngxmax, num_bins, 2, 2, ngxmax, ihymax+1,
     #     ngymax, num_bins, num_bins, ihymax+1, num_bins/),
     #    (/num_type_file-1,2/)) !values depending on *max may change below
                               !depending on the input file in use

      real, allocatable, dimension(:,:) ::   data_ar
      integer, allocatable, dimension(:,:) :: int_data_ar
      integer(kind=2), allocatable, dimension(:,:) :: out_int_data_ar
c Going to allocate the data array based on what type of file is being
c dealt with. Only one of data_ar or int_data_ar used in any particular case.

      real (kind=REAL_OUTPUT_KIND), allocatable, dimension(:,:) ::  
     #     out_data_ar
C  out_data_ar used to write the data to the hdf files.

      real(kind=REAL_OUTPUT_KIND), allocatable, dimension(:,:,:) ::  
     #	data_av_ar_hdf
      real, allocatable, dimension(:,:,:) ::  data_av_ar_unf

      real bin_bounds(2,2)  ! used only for type 3 - distribution funcs

C   data_av_ar used when the ave flag is set - its where intermediate
c   data which will be averaged down is set.
	
      integer,dimension(nsmax) :: attributes=0 ! only used with type 3      
c Going to allocate the data array based on what type of file is being
c dealt with.
c
      integer :: num_particles
      logical, optional :: average, drop_first_set
      integer, optional :: ave_over
      integer*4, optional :: sd_id_input

c      character*(maxchar), optional :: in_filename
C Ave will be set true in the case where this subroutine is supposed to
c do averaging of data.
  
c	write(*,*) "just entered unf ", filename, unf_unit, type,
c     #    data_name,  num_sets,it_0,it_stride, average

c  check to see if ihist and ihymax equal, if not change all_sizes values that
c depended on ihymax, because they are supposed to depend on ihist

      if (ihist .ne. ihymax) then
          all_sizes(1,2) = ihist+1
c	  all_sizes(5,2)=nsp*(ihist+1)
          all_sizes(5,2)=(ihist+1)
      endif



c and for ngxmax, ngymax, ngx, ngy
      if ((ngx .ne. ngxmax) .or. (ngy .ne. ngymax)) then
	  if (ngx .ne. ngxmax) then
	      all_sizes(2,1)=ngx
	      all_sizes(6,1)=ngx
	  else
	      all_sizes(2,2)=ngy
	  endif
      endif

 	   
  

C allocate the data array
      select case (type)
        ! every type except 3 and 6 has real data
        case(3,6)
          allocate (int_data_ar(all_sizes(type,1),all_sizes(type,2)))
          allocate(out_int_data_ar(all_sizes(type,1),all_sizes(type,2)))
        case default
          allocate (data_ar(all_sizes(type,1),all_sizes(type,2)))
          allocate (out_data_ar(all_sizes(type,1),all_sizes(type,2)))
          data_ar=12345.0
      end select

c seems confused, but..
c      allocate (data_ar(all_sizes(type,2),all_sizes(type,1)))


C  Make the definitions that will be the same for all file types.

C     Create and open the file and initiate the SD interface.
      sd_id = sfstart(filename(1:name_len), DFACC_WRITE)
      start(1) = 0
      start(2) = 0

C	Always start at 0 and stride 1.
C   stride can be used to make the data skip, 0 or 1 makes every space be used
      stride(1)=1
      stride(2)=1
      dimsizes(1) =all_sizes(type,1)
      dimsizes(2) =all_sizes(type,2)
      edges(1) = all_sizes(type,1)
      edges(2) = all_sizes(type,2)

      rank=2
c     data is either x vs y or contour type - both are rank 2.

c     The rest depends on data type - defined in setup.

      select case (type)  
        case(2)  !First, the contour case

c		write(*,*) "in case 2 "
C  allocate array for storing intermediate data if ave is set

	  if (present(average)) then

               if (present(sd_id_input)) then

			allocate (data_av_ar_hdf(ave_over, all_sizes(type,1),
     #		     		all_sizes(type,2)))
     		else
     			allocate (data_av_ar_unf(ave_over, all_sizes(type,1),
     #				all_sizes(type,2) ))
     		endif
		
		if (drop_first_set) then
			read(unf_unit, IOSTAT=ios_read_unf, ERR=1111,
     # 				END=1111) data_av_ar_unf(1,:,:)		        
		endif
C  The above line reads in the iteration 0 results for phi.  These results
C  are ignored, so that the averaging is like that in the main program.
C  Kinda ludgy adding an extra parameter (drop_first_set) for this little case,
C  but that's the way it goes I guess.

	  endif
	  
	  do i=1,num_sets
	    if (i .ne. num_sets) then
		it_num=it_0 +(i-1)*it_stride  ! figure out the iteration number
	    elseif (((unf_unit-unit_offset) .ne. 14) .and.((unf_unit
     #	    	- unit_offset) .ne. 0) ) then 
		it_num=nt  ! always output on last iteration
	    else
		it_num=nt-it_0 ! special case for phiav or phiave because
		    ! average is for middle of range - not real sure if
		    ! this will always work.
	    endif

c	    write(*,*) nt,it_num 
            call int2char_p(it_num,it_num_ch)
	    set_label=data_name(1:dat_name_len)//": It. "//
     #         it_num_ch
c            write(*,*) it_num,it_num_ch,set_label
C  First, read the unformatted data into an array.
c
C  For speed would read entire file in at once.  Reading in set by set saves
C  memory.

C  If we're avergaing down the data, have to manipulate it differently.

	    if (present(average)) then

c	   	if (present(in_filename)) then   ! reading from an hdf file
               if (present(sd_id_input)) then


			do j=1,ave_over
			    call read_hdf_dataset( (i-1)
     # 				*ave_over + j, all_sizes(2, :),
     #				data_av_ar_hdf(j,:,:), sd_id_input)


C           Average values (using spiffy matrix SUM operator)
C	    Its kind of kludgy to change the kind of data_ar now and later,
C	    but its not worth changing at this point.

		            data_ar = REAL((SUM (data_av_ar_hdf, 1)
     #				/ave_over), DEFAULT_KIND)
				 !sum array over first dim

     			    
     			enddo
			
	     	else   !reading from unformatted file

			do j=1,ave_over
	    		    read(unf_unit, IOSTAT=ios_read_unf,
     #			    	ERR=1111,END=1111)data_av_ar_unf(j,:,:)

C           Average values (using spiffy matrix SUM operator)
		            data_ar = SUM(data_av_ar_unf, 1)/ave_over
				 !sum array over first dim

 
     			enddo
			
		endif  ! if present
		     
            else    ! if not averaging, simply read in a set.
    	    	read(unf_unit,IOSTAT=ios_read_unf, ERR=1111, 
     #			END=1111) data_ar
	    endif

c	    write(*,*) "about to create a dataset"

	    out_data_ar = REAL(data_ar, REAL_OUTPUT_KIND)
C  Converts data from the default real type, to the chosen real type.  This
C  uses a little CPU time, but saves alot of disk space.

            sds_id=sfcreate(sd_id, set_label, REAL_HDF_KIND, rank,
     #	    	dimsizes)
	    status = sfwdata(sds_id, start, stride,edges, out_data_ar)
	    status = sfendacc(sds_id)
          enddo

        case(0:1) ! 
	  !write(*,*) 'in case 1', num_sets
          do i=1,num_sets
	    if (i .ne. num_sets) then
                it_num=it_0 +(i-1)*it_stride  ! figure out the iteration number
            else 
                it_num=nt  ! always output on last iteration 
	    endif

            call int2char_p(it_num,it_num_ch)
            set_label=data_name(1:dat_name_len)//": It. "//
     #         it_num_ch
            read(unf_unit, IOSTAT=ios_read_unf, ERR=1111, END=1111)
     #	    	data_ar(1,1:all_sizes(type,2)), data_ar(2,1:
     #		all_sizes(type,2))
c This separation of the read is necesary because that's the way the data is
c written - 1 column at a time. Luckily, this shouldn't be much of a performance
c problem since columns of data are stored contigously anyway.

	    out_data_ar = REAL(data_ar, REAL_OUTPUT_KIND)
            sds_id=sfcreate(sd_id,set_label,REAL_HDF_KIND,rank,dimsizes)
            status = sfwdata(sds_id, start, stride,edges, out_data_ar)
            status = sfendacc(sds_id)
          enddo

       case(3,6)
c          write(*,*) "in case 3"
          do i=1,num_sets
            if (i .ne. num_sets) then
                it_num=it_0 +(i-1)*it_stride  ! figure out the iteration number
            else
                it_num=nt  ! always output on last iteration 
            endif

            call int2char_p(it_num,it_num_ch)   ! make iteration number into
	    call getlen(it_num_ch,it_num_len)   ! a character variable
            !it_num_len = len(it_num_ch)

            do jj=1,nsp  ! species dependant case
                call int2char_p(jj,spec_char)  ! get char for species number
                set_label=data_name(1:dat_name_len)//": It. "//
     #           	it_num_ch(1:it_num_len)//": Sp. "//spec_char

                read(unf_unit, IOSTAT=ios_read_unf, ERR=1111, END=1111)
     #	    		bin_bounds

            	read(unf_unit, IOSTAT=ios_read_unf, ERR=1111, END=1111)
     #	    		int_data_ar


                sds_id=sfcreate(sd_id,set_label,DFNT_INT16,rank,
     #			dimsizes)
	    	! Should do some type of error catching here.  If one of these 
	    	! integers ends up being greater than 32k,there'll be problems.
	    	out_int_data_ar = int(int_data_ar, 2) 
		!reduce the kind to save space

c Set the maximum and minimum as a separate attribute for each.
c The three is the size of the vmin array.
            	status=sfsnatt( sds_id, 'Min and max ',REAL_HDF_KIND,
     #	    		4, REAL(bin_bounds, REAL_OUTPUT_KIND)) 

            	status = sfwdata(sds_id, start, stride, edges,
     #			out_int_data_ar)
	    	!status = sfwdata(sds_id, start, stride, edges, int_data_ar)
            	status = sfendacc(sds_id)
	    enddo
          enddo

        case(4:5)
c	  write(*,*) "in case 4 or 5"
          do i=1,num_sets
            if (i .ne. num_sets) then
                  it_num=it_0 +(i-1)*it_stride  !figure out the iteration number
            else
                it_num=nt  ! always output on last iteration 
            endif
            call int2char_p(it_num,it_num_ch)
	    call getlen(it_num_ch,it_num_len)
	    !it_num_len = len(it_num_ch)

C		write(*,*) num_sets,it_0,it_stride,it_num,'...'

            do jj=1,nsp  ! species dependant case
              call int2char_p(jj,spec_char)  ! get char for species number
              set_label=data_name(1:dat_name_len)//": It. "//
     #           it_num_ch(1:it_num_len)//": Sp. "//spec_char
              read(unf_unit,IOSTAT=ios_read_unf, ERR=1111, END=1111)
     #	      	 data_ar(1,1:all_sizes(type,2)), data_ar(2,1:
     #		 all_sizes(type,2))
c This separation of the read is necessary because that's the way the data is
c written - 1 column at a time. Luckily, this shouldn't be much of a performance
c problem since columns of data are stored contigously anyway.

	      out_data_ar = REAL(data_ar, REAL_OUTPUT_KIND)

              sds_id=sfcreate(sd_id,set_label,REAL_HDF_KIND,rank,
     #           dimsizes)
              status = sfwdata(sds_id, start, stride,edges, out_data_ar)
              status = sfendacc(sds_id)
         
	    enddo

	  enddo

	case default
          write(*,*) 'Found a bad choice in unf_to_hdf - oops!', type

      end select

C The next line is jumped to if ther is an error in reading one of the unf 
C files.
 1111 if (ios_read_unf .ne. 0) then   ! Check the error condition
 		
		write(*,*) "There was an error opening the binary ",
     #			"file associated with ", filename(1:name_len),
     #			" for reading. Hopefully you can figure out ", 
     #			"the cause. The error occurred in the ",
     # 			" unf_to_hdf routine. Have a wonderful day.",
     #			i, jj, ios_read_unf
		ios_read_unf = 0  ! reset the value, so error message doesn't
					! pop up for the rest of files.

      endif

C Stuff below here should be the same for all cases.

C     Terminate access to the array.
 
C     Terminate access to the SD interface and close the file.
      status = sfend(sd_id)

C     Close the unformatted data file.
      close (unf_unit)

c Free up the space for the arrays
      select case (type)  

	case(3,6)
               deallocate (int_data_ar, out_int_data_ar)

	case default
                deallocate (data_ar, out_data_ar)

      end select


      if (present(average)) then
c                if (present(in_filename)) then
                if (present(sd_id_input)) then
                        deallocate (data_av_ar_hdf)
                else
                        deallocate (data_av_ar_unf)
                endif

      endif

c	write(*,*) "leaving unf"
      end  subroutine unf_to_hdf



      end module unf_to_hdf_mod
 

