output_collapsed_vti Subroutine

subroutine output_collapsed_vti(dir,normconv)
use mod_global_parameters
use mod_calculate_xw
integer, intent(in)                               :: dir
double precision,dimension(0:nw+nwauxio),intent(in):: normconv 
character(len=1024) :: filename, outfilehead, line
logical             :: fileopen
integer                                           :: iw
character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
integer, dimension(ndim)                          :: myshape

integer                                           :: ix1,ix2
double precision                                  :: dxdim1,dxdim2, xdimmin1,&
   xdimmin2,xdimmax1,xdimmax2

double precision                                  :: origin(1:3), spacing(1:3)
integer                                           :: wholeExtent(1:6),&
    size_single, length, size_length
integer*8                                         :: offset
character                                         :: buf
!-----------------------------------------------------------------------------
if (mype==0) then
offset=0
size_single=4
size_length=4
! Get coordinates:

select case(dir)
case (1)
   dxdim1 = dx(2,1)*2.0d0**(1-collapseLevel)
   dxdim2 = dx(3,1)*2.0d0**(1-collapseLevel)
   xdimmin1=xprobmin2;xdimmax1=xprobmax2; 
   xdimmin2=xprobmin3;xdimmax2=xprobmax3; 
case (2)
   dxdim1 = dx(1,1)*2.0d0**(1-collapseLevel)
   dxdim2 = dx(3,1)*2.0d0**(1-collapseLevel)
   xdimmin1=xprobmin1;xdimmax1=xprobmax1; 
   xdimmin2=xprobmin3;xdimmax2=xprobmax3; 
case (3)
   dxdim1 = dx(1,1)*2.0d0**(1-collapseLevel)
   dxdim2 = dx(2,1)*2.0d0**(1-collapseLevel)
   xdimmin1=xprobmin1;xdimmax1=xprobmax1; 
   xdimmin2=xprobmin2;xdimmax2=xprobmax2; 
case default
   dxdim1 = 1
   dxdim2 = 1
   xdimmin1=xprobmin2;xdimmax1=xprobmax2; 
   xdimmin2=xprobmin3;xdimmax2=xprobmax3; 
   call mpistop("slice direction not clear in output_collapsed_vti")
end select



origin=0
spacing=zero
wholeExtent=0
myshape = shape(collapsedData)


length = myshape(1)*myshape(2)
length = length*size_single
wholeExtent(1*2)=myshape(1); wholeExtent(2*2)=myshape(2);
spacing(1) = dxdim1; spacing(2) = dxdim2;
origin(1) = xdimmin1; origin(2) = xdimmin2;


 inquire(unitcollapse,opened=fileopen)
 if(.not.fileopen)then
      ! generate filename: 
    write(filename,"(a,i1.1,a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,&
       '_l',collapseLevel,'_n',collapsenext,'.vti'
      open(unitcollapse,file=filename,status='unknown',form='formatted')
 end if
! get the header: 
call getheadernames(wnamei,xandwnamei,outfilehead)

! generate xml header
write(unitcollapse,'(a)')'<?xml version="1.0"?>'
write(unitcollapse,'(a)',advance='no') '<VTKFile type="ImageData"'
if(type_endian==1)then
   write(unitcollapse,'(a)')' version="0.1" byte_order="LittleEndian">'
else
  write(unitcollapse,'(a)')' version="0.1" byte_order="BigEndian">'
endif
! following corresponds to uniform cartesian grid
write(unitcollapse,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')&
   '  <ImageData Origin="',origin,'" WholeExtent="',wholeExtent,'" Spacing="',&
   spacing,'">'
write(unitcollapse,'(a)')'<FieldData>'
write(unitcollapse,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
   'NumberOfTuples="1" format="ascii">'
write(unitcollapse,*) real(global_time*time_convert_factor)
write(unitcollapse,'(a)')'</DataArray>'
write(unitcollapse,'(a)')'</FieldData>'

! we write one VTK PIECE
write(unitcollapse,'(a,6(i10),a)') '<Piece Extent="',wholeExtent,'">'
write(unitcollapse,'(a)')'<CellData>'

do iw=1,nw+nwauxio
  if(iw<=nw) then 
    if(.not.w_write(iw)) cycle
  endif
  write(unitcollapse,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
     trim(wnamei(iw)),'" format="appended" offset="',offset,'"/>'
  offset = offset + length + size_length
enddo

write(unitcollapse,'(a)')'</CellData>'
write(unitcollapse,'(a)')'</Piece>'
write(unitcollapse,'(a)')'</ImageData>'
write(unitcollapse,'(a)')'<AppendedData encoding="raw">'
close(unitcollapse)

open(unitcollapse,file=filename,access='stream',form='unformatted',&
   position='append')
buf='_'
write(unitcollapse) TRIM(buf)

do iw=1,nw+nwauxio
  if(iw<=nw) then 
    if(.not.w_write(iw)) cycle
  endif
  write(unitcollapse) length
  
  
  
  write(unitcollapse) ((real(collapsedData(ix1,ix2,iw)*normconv(iw)),ix1=1,&
     myshape(1)),ix2=1,myshape(2))
 
enddo

close(unitcollapse)
open(unitcollapse,file=filename,status='unknown',form='formatted',&
   position='append')
write(unitcollapse,'(a)')'</AppendedData>'
write(unitcollapse,'(a)')'</VTKFile>'
close(unitcollapse)
 
end if

end subroutine output_collapsed_vti