subroutine put_slice(dir,xslice)
use mod_forest, only: Morton_sub_start, Morton_sub_stop
use mod_global_parameters
use mod_geometry, only: s_of_r
! Writes a D-1 slice
! For ONED simulations, output will be appended to one csv-file per slice
! slices are sensitive to the saveprim switch and
! can contain auxiliary io variables (nwauxio)
! Thus csv or vtu(CC)-files with primitive variables are obtained.
! when slice_type='dat', we save a D-1 dat file for potential restarts
! by Oliver Porth
! 22.Nov 2011
integer, intent(in) :: dir
double precision, intent(in) :: xslice
! .. local ..
integer :: Njgrid, jgrid
integer, dimension(ndim-1) :: ixsubGlo, ixsubGhi
integer, dimension(ndim-1) :: ixsubMlo, ixsubMhi
integer :: size_subblock_io, nx1,nx2,nx3, slice_fh, nwexpand
integer :: type_subblock_io, type_subblockC_io, type_subblock_x_io,&
type_subblockC_x_io
integer, dimension(ndim) :: sizes, subsizes, start
double precision,dimension(0:nw+nwauxio) :: normconv
double precision, dimension(ndim-1) :: xprobminsub, xprobmaxsub
integer, dimension(ndim-1) :: block_nx_sub, domain_nx_sub
logical, dimension(ndim-1) :: periodBsub
! Preamble:
nx1=ixMhi1-ixMlo1+1;nx2=ixMhi2-ixMlo2+1;nx3=ixMhi3-ixMlo3+1;
slice_fh=unitslice
if (ndim==1) then
nwexpand = nwauxio
else
if(slice_type/='dat')then
nwexpand = nwauxio
else
nwexpand = 0
endif
end if
! Do a last consistency check. xslice is a physical coordinate, so under
! LOG_RADIUS the radial one has to be compared in the logical coordinate
! that xprobmin1/xprobmax1 are stored in.
select case(dir)
case(1)
if(xslice<xprobmin1.or.xslice>xprobmax1) call &
mpistop("slice out of bounds")
case(2)
if(xslice<xprobmin2.or.xslice>xprobmax2) call &
mpistop("slice out of bounds")
case(3)
if(xslice<xprobmin3.or.xslice>xprobmax3) call &
mpistop("slice out of bounds")
end select
! Traverse the forest and fill nodes:
call select_slice(dir,xslice,.false.,slice_fh,normconv)
! Create the MPI-datatype and select indices:
select case(dir)
case (1)
sizes(1) = ixGhi2; sizes(2) = ixGhi3;
ixsubGlo(1) = ixGlo2; ixsubGlo(2) = ixGlo3;
ixsubGhi(1) = ixGhi2; ixsubGhi(2) = ixGhi3;
subsizes(1)=nx2;subsizes(2)=nx3;
start(1)=ixMlo2-1;start(2)=ixMlo3-1;
size_subblock_io=nx2*nx3*(nw+nwexpand)*size_double
xprobminsub(1)=xprobmin2; xprobmaxsub(1)=xprobmax2;
xprobminsub(2)=xprobmin3; xprobmaxsub(2)=xprobmax3;
domain_nx_sub(1)=domain_nx2; domain_nx_sub(2)=domain_nx3;
block_nx_sub(1)=block_nx2; block_nx_sub(2)=block_nx3;
periodBsub(1)=periodB(2); periodBsub(2)=periodB(3);
case (2)
sizes(1) = ixGhi1; sizes(2) = ixGhi3;
ixsubGlo(1) = ixGlo1; ixsubGlo(2) = ixGlo3;
ixsubGhi(1) = ixGhi1; ixsubGhi(2) = ixGhi3;
subsizes(1)=nx1;subsizes(2)=nx3;
start(1)=ixMlo1-1;start(2)=ixMlo3-1;
size_subblock_io=nx1*nx3*(nw+nwexpand)*size_double
xprobminsub(1)=xprobmin1; xprobmaxsub(1)=xprobmax1;
xprobminsub(2)=xprobmin3; xprobmaxsub(2)=xprobmax3;
domain_nx_sub(1)=domain_nx1; domain_nx_sub(2)=domain_nx3;
block_nx_sub(1)=block_nx1; block_nx_sub(2)=block_nx3;
periodBsub(1)=periodB(1); periodBsub(2)=periodB(3);
case (3)
ixsubGlo(1) = ixGlo1; ixsubGlo(2) = ixGlo2;
ixsubGhi(1) = ixGhi1; ixsubGhi(2) = ixGhi2;
sizes(1) = ixGhi1; sizes(2) = ixGhi2;
subsizes(1)=nx1;subsizes(2)=nx2;
start(1)=ixMlo1-1;start(2)=ixMlo2-1;
size_subblock_io=nx1*nx2*(nw+nwexpand)*size_double
xprobminsub(1)=xprobmin1; xprobmaxsub(1)=xprobmax1;
xprobminsub(2)=xprobmin2; xprobmaxsub(2)=xprobmax2;
domain_nx_sub(1)=domain_nx1; domain_nx_sub(2)=domain_nx2;
block_nx_sub(1)=block_nx1; block_nx_sub(2)=block_nx2;
periodBsub(1)=periodB(1); periodBsub(2)=periodB(2);
case default
call mpistop("slice direction not clear in put_slice")
end select
ixsubMlo(2-1) = ixsubGlo(2-1)+nghostcells
ixsubMlo(3-1) = ixsubGlo(3-1)+nghostcells;
ixsubMhi(2-1) = ixsubGhi(2-1)-nghostcells
ixsubMhi(3-1) = ixsubGhi(3-1)-nghostcells;
sizes(ndim)=nw+nwexpand
subsizes(ndim)=nw+nwexpand
start(ndim)=0
! Types for center variables:
call MPI_TYPE_CREATE_SUBARRAY(ndim,sizes,subsizes,start, MPI_ORDER_FORTRAN,&
MPI_DOUBLE_PRECISION, type_subblock_io,ierrmpi)
call MPI_TYPE_COMMIT(type_subblock_io,ierrmpi)
sizes(ndim)=3
subsizes(ndim)=3
start(ndim)=0
call MPI_TYPE_CREATE_SUBARRAY(ndim,sizes,subsizes,start, MPI_ORDER_FORTRAN,&
MPI_DOUBLE_PRECISION, type_subblock_x_io,ierrmpi)
call MPI_TYPE_COMMIT(type_subblock_x_io,ierrmpi)
! Types for corner variables:
subsizes(1:ndim-1) = subsizes(1:ndim-1) + 1
start(1:ndim-1) = start(1:ndim-1) - 1
sizes(ndim)=nw+nwexpand
subsizes(ndim)=nw+nwexpand
start(ndim)=0
call MPI_TYPE_CREATE_SUBARRAY(ndim,sizes,subsizes,start, MPI_ORDER_FORTRAN,&
MPI_DOUBLE_PRECISION, type_subblockC_io,ierrmpi)
call MPI_TYPE_COMMIT(type_subblockC_io,ierrmpi)
sizes(ndim)=3
subsizes(ndim)=3
start(ndim)=0
call MPI_TYPE_CREATE_SUBARRAY(ndim,sizes,subsizes,start, MPI_ORDER_FORTRAN,&
MPI_DOUBLE_PRECISION, type_subblockC_x_io,ierrmpi)
call MPI_TYPE_COMMIT(type_subblockC_x_io,ierrmpi)
! local number of sub-grids:
Njgrid=Morton_sub_stop(mype)-Morton_sub_start(mype)+1
! Now output using various schemes:
if (ndim==1) then
call put_slice_zerod
else
select case(slice_type)
case ('csv')
call put_slice_csv
case ('dat')
call put_slice_dat
case ('vtu', 'vtuCC')
call put_slice_vtu
end select
end if
do jgrid=1,Njgrid
call dealloc_subnode(jgrid)
end do
call MPI_TYPE_FREE(type_subblock_io,ierrmpi)
call MPI_TYPE_FREE(type_subblock_x_io,ierrmpi)
call MPI_TYPE_FREE(type_subblockC_io,ierrmpi)
call MPI_TYPE_FREE(type_subblockC_x_io,ierrmpi)
if (allocated(sfc_sub)) deallocate(sfc_sub)
contains
subroutine put_slice_vtu
use mod_calculate_xw
character(len=1024) :: filename, xlabel
character(len=79) :: xxlabel
logical :: fileopen
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
integer :: status(MPI_STATUS_SIZE), ipe
if (mype==0) then
inquire(slice_fh,opened=fileopen)
if(.not.fileopen)then
! generate filename:
write(xlabel,"(D9.2)")xslice
xxlabel=trim(xlabel)
if(xslice>=zero)then
write(xxlabel(1:1),"(a)") "+"
endif
write(filename,"(a,i1.1,a,i4.4,a)") TRIM(base_filename)//'_d',dir,&
'_x'//trim(xxlabel)//'_n',slicenext,'.vtu'
open(slice_fh,file=filename,status='unknown',form='formatted')
end if
! get and write the header:
call getheadernames(wnamei,xandwnamei,outfilehead)
! generate xml header
write(slice_fh,'(a)')'<?xml version="1.0"?>'
write(slice_fh,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
if(type_endian==1)then
write(slice_fh,'(a)')' version="0.1" byte_order="LittleEndian">'
else
write(slice_fh,'(a)')' version="0.1" byte_order="BigEndian">'
endif
write(slice_fh,'(a)')' <UnstructuredGrid>'
write(slice_fh,'(a)')'<FieldData>'
write(slice_fh,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
'NumberOfTuples="1" format="ascii">'
write(slice_fh,*) real(global_time*time_convert_factor)
write(slice_fh,'(a)')'</DataArray>'
write(slice_fh,'(a)')'</FieldData>'
! write to file:
do jgrid=1, Njgrid
call write_slice_vtk(jgrid,slice_fh,wnamei)
end do
! create a recv buffer using allocate, will be deallocated at the end of the routine:
call alloc_subnode(Njgrid+1,dir,nwauxio)
end if
! Also communicate the normconv array since processor zero might not have it yet:
if (npe>1) then
do ipe=1,npe-1
do jgrid=1,Morton_sub_stop(ipe)-Morton_sub_start(ipe)+1
itag=Morton_sub_start(ipe)+jgrid-1
itag=itag*5
if (ipe == mype ) then
call mpi_send_wrapper(ps_sub(jgrid)%x,1,type_subblock_x_io,0,itag,&
icomm,ierrmpi)
call mpi_send_wrapper(ps_sub(jgrid)%w,1,type_subblock_io,0,itag+1,&
icomm,ierrmpi)
call mpi_send_wrapper(ps_sub(jgrid)%xC,1,type_subblockC_x_io,0,&
itag+2,icomm,ierrmpi)
call mpi_send_wrapper(ps_sub(jgrid)%wC,1,type_subblockC_io,0,itag+3,&
icomm,ierrmpi)
call mpi_send_wrapper(normconv,nw+nwauxio+1,MPI_DOUBLE_PRECISION,0,&
itag+4,icomm,ierrmpi)
end if
if (mype == 0) then
call mpi_recv_wrapper(ps_sub(Njgrid+1)%x,1,type_subblock_x_io,ipe,&
itag,icomm,status,ierrmpi)
call mpi_recv_wrapper(ps_sub(Njgrid+1)%w,1,type_subblock_io,ipe,&
itag+1,icomm,status,ierrmpi)
call mpi_recv_wrapper(ps_sub(Njgrid+1)%xC,1,type_subblockC_x_io,ipe,&
itag+2,icomm,status,ierrmpi)
call mpi_recv_wrapper(ps_sub(Njgrid+1)%wC,1,type_subblockC_io,ipe,&
itag+3,icomm,status,ierrmpi)
call mpi_recv_wrapper(normconv,nw+nwauxio+1,MPI_DOUBLE_PRECISION,ipe,&
itag+4,icomm,status,ierrmpi)
call write_slice_vtk(Njgrid+1,slice_fh,wnamei)
end if
end do
end do
endif
if (mype==0) then
write(slice_fh,'(a)')'</UnstructuredGrid>'
write(slice_fh,'(a)')'</VTKFile>'
close(slice_fh)
call dealloc_subnode(Njgrid+1)
end if
end subroutine put_slice_vtu
subroutine write_slice_vtk(jgrid,slice_fh,wnamei)
! this only works for 2D and 3D, 1D reduction (line to point) not allowed
integer, intent(in) :: jgrid, slice_fh
character(len=name_len), intent(in) :: wnamei(1:nw+nwauxio)
! This remainder part only for more than 1D, but nesting with NOONED gives problems
! .. local ..
integer :: ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,&
ixCmax3, ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3, nc,&
np, iw
integer :: nx1,nx2, nxC1,nxC2, icell, ix1,ix2
double precision :: x_VTK(1:3)
integer :: VTK_type
double precision, parameter :: minvalue = 1.0d-99, maxvalue = 1.0d+99
ixCCmin1 = ixsubMlo(1);ixCCmin2 = ixsubMlo(2);
ixCCmax1 = ixsubMhi(1);ixCCmax2 = ixsubMhi(2);
ixCmin1 = ixsubMlo(1)-1;ixCmin2 = ixsubMlo(2)-1;
ixCmax1 = ixsubMhi(1);ixCmax2 = ixsubMhi(2);
nx1=ixCCmax1-ixCCmin1+1;nx2=ixCCmax2-ixCCmin2+1;
nxC1=nx1+1;nxC2=nx2+1;
nc=nx1*nx2 ! Number of cells per subgrid
np=nxC1*nxC2 ! Number of corner points per subgrid
! we write out every grid as one VTK PIECE
write(slice_fh,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
'" NumberOfCells="',nc,'">'
!==============================
! celldata or pointdata?
!==============================
select case(slice_type)
case('vtuCC') ! celldata
write(slice_fh,'(a)')'<CellData>'
do iw=1,nw+nwauxio
if(iw<=nw) then
if(.not.w_write(iw)) cycle
endif
write(slice_fh,'(a,a,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)),'" format="ascii">'
write(slice_fh,'(200(1pe14.6))') &
((roundoff_minmax(ps_sub(jgrid)%w(ix1,ix2,iw)*normconv(iw),&
minvalue,maxvalue),ix1=ixCCmin1,ixCCmax1),ix2=ixCCmin2,&
ixCCmax2)
write(slice_fh,'(a)')'</DataArray>'
enddo
write(slice_fh,'(a)')'</CellData>'
case('vtu') ! pointdata
write(slice_fh,'(a)')'<PointData>'
do iw=1,nw+nwauxio
if(iw<=nw) then
if(.not.w_write(iw)) cycle
endif
write(slice_fh,'(a,a,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)),'" format="ascii">'
write(slice_fh,'(200(1pe14.6))') &
((roundoff_minmax(ps_sub(jgrid)%wC(ix1,ix2,iw)*normconv(iw),&
minvalue,maxvalue),ix1=ixCmin1,ixCmax1),ix2=ixCmin2,ixCmax2)
write(slice_fh,'(a)')'</DataArray>'
enddo
write(slice_fh,'(a)')'</PointData>'
end select
!==============================
! Done: celldata or pointdata?
!==============================
!==============================
! output Cornerpoints
!==============================
write(slice_fh,'(a)')'<Points>'
write(slice_fh,'(a)'&
)'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
!write cell corner coordinates in a backward dimensional loop, always 3D output
do ix2=ixCmin2,ixCmax2
do ix1=ixCmin1,ixCmax1
x_VTK(1:3)=zero;
x_VTK(1:ndim)=ps_sub(jgrid)%xC(ix1,ix2,1:ndim)*normconv(0);
write(slice_fh,'(3(1pe14.6))') x_VTK
end do
end do
write(slice_fh,'(a)')'</DataArray>'
write(slice_fh,'(a)')'</Points>'
!==============================
! Done: output Cornerpoints
!==============================
!==============================
! cell Metainformation
!==============================
write(slice_fh,'(a)')'<Cells>'
! connectivity part
write(slice_fh,'(a)'&
)'<DataArray type="Int32" Name="connectivity" format="ascii">'
do ix2=1,nx2
do ix1=1,nx1
write(slice_fh,'(4(i7,1x))')(ix2-1)*nxC1+ix1-1, (ix2-1)*nxC1+ix1,&
ix2*nxC1+ix1-1,ix2*nxC1+ix1
end do
end do
write(slice_fh,'(a)')'</DataArray>'
! offsets data array
write(slice_fh,'(a)'&
)'<DataArray type="Int32" Name="offsets" format="ascii">'
do icell=1,nc
write(slice_fh,'(i7)') icell*(2**(3-1))
end do
write(slice_fh,'(a)')'</DataArray>'
! VTK cell type data array
write(slice_fh,'(a)'&
)'<DataArray type="Int32" Name="types" format="ascii">'
! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
VTK_type=8
do icell=1,nc
write(slice_fh,'(i2)') VTK_type
enddo
write(slice_fh,'(a)')'</DataArray>'
write(slice_fh,'(a)')'</Cells>'
!==============================
! Done: cell Metainformation
!==============================
write(slice_fh,'(a)')'</Piece>'
end subroutine write_slice_vtk
subroutine put_slice_csv
use mod_calculate_xw
character(len=1024) :: filename, xlabel
character(len=79) :: xxlabel
logical :: fileopen
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
integer :: iw, ipe, itag
character(len=1024) :: line
integer :: status(MPI_STATUS_SIZE)
if (mype==0) then
inquire(slice_fh,opened=fileopen)
if(.not.fileopen)then
! generate filename:
write(xlabel,"(D9.2)")xslice
xxlabel=trim(xlabel)
if(xslice>=zero)then
write(xxlabel(1:1),"(a)") "+"
endif
write(filename,"(a,i1.1,a,i4.4,a)") TRIM(base_filename)//'_d',dir,&
'_x'//trim(xxlabel)//'_n',slicenext,'.csv'
open(slice_fh,file=filename,status='unknown',form='formatted')
end if
! get and write the header:
call getheadernames(wnamei,xandwnamei,outfilehead)
line=''
do iw=1,ndim+nw+nwauxio-1
line = trim(line)//trim(xandwnamei(iw))//', '
end do
line = trim(line)//trim(xandwnamei(ndim+nw+nwauxio))
write(slice_fh,'(a)')trim(line)
! create a recv buffer using allocate, will be deallocated at the end of the routine:
call alloc_subnode(Njgrid+1,dir,nwauxio)
! write to file:
do jgrid=1, Njgrid
call put_slice_line(jgrid,slice_fh)
end do
end if
! Also communicate the normconv array since processor zero might not have it yet:
if (npe>1) then
do ipe=1,npe-1
do jgrid=1,Morton_sub_stop(ipe)-Morton_sub_start(ipe)+1
itag=Morton_sub_start(ipe)+jgrid-1
if (ipe == mype ) then
call mpi_send_wrapper(ps_sub(jgrid)%x,1,type_subblock_x_io,0,itag,&
icomm,ierrmpi)
call mpi_send_wrapper(ps_sub(jgrid)%w,1,type_subblock_io,0,itag,&
icomm,ierrmpi)
call mpi_send_wrapper(normconv,nw+nwauxio+1,MPI_DOUBLE_PRECISION,0,&
itag,icomm,ierrmpi)
end if
if (mype == 0) then
call mpi_recv_wrapper(ps_sub(Njgrid+1)%x,1,type_subblock_x_io,ipe,&
itag,icomm,status,ierrmpi)
call mpi_recv_wrapper(ps_sub(Njgrid+1)%w,1,type_subblock_io,ipe,itag,&
icomm,status,ierrmpi)
call mpi_recv_wrapper(normconv,nw+nwauxio+1,MPI_DOUBLE_PRECISION,ipe,&
itag,icomm,status,ierrmpi)
call put_slice_line(Njgrid+1,slice_fh)
end if
end do
end do
endif
if (mype==0) then
close(slice_fh)
call dealloc_subnode(Njgrid+1)
end if
end subroutine put_slice_csv
subroutine put_slice_line(jout,file_handle)
integer, intent(in) :: jout, file_handle
! .. local ..
character(len=1024) ::line, data
integer :: ix1,ix2,ix3,idir,iw
double precision, parameter :: minvalue = 1.0d-99, maxvalue = 1.0d+99
do ix2=ixsubMlo(2),ixsubMhi(2)
do ix1=ixsubMlo(1),ixsubMhi(1)
! Format the line:
line = ''
do idir=1,ndim
write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(ix1,ix2,&
idir),minvalue,maxvalue)
line = trim(line)//trim(data)//', '
end do
do iw = 1,nw+nwauxio-1
write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,ix2,&
iw)*normconv(iw),minvalue,maxvalue)
line = trim(line)//trim(data)//', '
end do
write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,ix2,&
nw+nwauxio)*normconv(nw+nwauxio),minvalue,maxvalue)
line = trim(line)//trim(data)
write(file_handle,'(a)')trim(line)
end do
end do
end subroutine put_slice_line
subroutine put_slice_dat
use mod_forest, only: tree_node, igrid_to_node
use mod_input_output_helper, only: version_number, save_now
use mod_physics, only: physics_type, phys_write_info
integer, dimension(max_blocks*2) :: iorequest
integer, dimension(MPI_STATUS_SIZE,max_blocks*2) :: iostatus
integer(kind=MPI_OFFSET_KIND) :: offset, offset1
integer :: nsubleafs
character(len=1024) :: filename, xlabel
character(len=79) :: xxlabel
integer :: amode, status(MPI_STATUS_SIZE), iwrite
character(len=name_len) :: dname
integer :: nsubparents, iw, ileaf, inode
integer :: n_ghost_sub(2*(ndim-1))
integer(kind=MPI_OFFSET_KIND) :: block_offset, tree_offset
integer, allocatable :: block_ig_sub(:,:)
integer, allocatable :: block_lvl_sub(:)
integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset_sub(:)
logical, allocatable :: isleaf_sub(:)
type(tree_node), pointer :: pnode
nsubleafs=Morton_sub_stop(npe-1)
if (mype==0) then
if (nsubleafs /= sfc_sub(0,1)) error stop "nsubleafs /= sfc_sub(0,1)"
nsubparents = sfc_sub(0,2)
allocate(block_ig_sub(ndim-1,nsubleafs))
allocate(block_lvl_sub(nsubleafs))
allocate(block_offset_sub(nsubleafs))
allocate(isleaf_sub(nsubleafs+nsubparents))
block_offset = int(0,kind=MPI_OFFSET_KIND)
ileaf = 0
do inode=1,nsubleafs+nsubparents
if (sfc_sub(inode,1) == 0) then
isleaf_sub(inode) = .false.
else
isleaf_sub(inode) = .true.
pnode => igrid_to_node(sfc_sub(inode,1),sfc_sub(inode,2))%node
ileaf = ileaf + 1
select case (dir)
case (1)
block_ig_sub(:,ileaf) = [ pnode%ig2, pnode%ig3 ]
case (2)
block_ig_sub(:,ileaf) = [ pnode%ig1, pnode%ig3 ]
case (3)
block_ig_sub(:,ileaf) = [ pnode%ig1, pnode%ig2 ]
end select
block_lvl_sub(ileaf) = pnode%level
block_offset_sub(ileaf) = block_offset
block_offset = block_offset + int(2*(ndim-1)*size_int,kind=MPI_OFFSET_KIND)
block_offset = block_offset + int(product(block_nx_sub)*nw*size_double,&
kind=MPI_OFFSET_KIND)
end if
end do
end if
! generate filename
write(xlabel,"(D9.2)")xslice
xxlabel=trim(xlabel)
if(xslice>=zero)then
write(xxlabel(1:1),"(a)") "+"
endif
write(filename,"(a,i1.1,a,i4.4,a)") TRIM(base_filename)//'_d',dir,&
'_x'//trim(xxlabel)//'_n',slicenext,'.dat'
if(mype==0) then
open(unit=slice_fh,file=filename,status='replace')
close(unit=slice_fh)
end if
if (mype==0) then
amode=ior(MPI_MODE_APPEND,MPI_MODE_WRONLY)
call MPI_FILE_OPEN(MPI_COMM_SELF,filename,amode,MPI_INFO_NULL,&
slice_fh,ierrmpi)
call MPI_FILE_WRITE(slice_fh, version_number, 1, MPI_INTEGER, status,&
ierrmpi)
call MPI_FILE_WRITE(slice_fh, 0, 1, MPI_INTEGER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, 0, 1, MPI_INTEGER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, nw, 1, MPI_INTEGER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, ndir, 1, MPI_INTEGER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, ndim-1, 1, MPI_INTEGER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, levmax_sub, 1, MPI_INTEGER, status,&
ierrmpi)
call MPI_FILE_WRITE(slice_fh, nsubleafs, 1, MPI_INTEGER, status,&
ierrmpi)
call MPI_FILE_WRITE(slice_fh, nsubparents, 1, MPI_INTEGER, status,&
ierrmpi)
call MPI_FILE_WRITE(slice_fh, it, 1, MPI_INTEGER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, global_time, 1, MPI_DOUBLE_PRECISION,&
status,ierrmpi)
call MPI_FILE_WRITE(slice_fh, xprobminsub, ndim-1,&
MPI_DOUBLE_PRECISION, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, xprobmaxsub, ndim-1,&
MPI_DOUBLE_PRECISION, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, domain_nx_sub, ndim-1, MPI_INTEGER,&
status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, block_nx_sub, ndim-1, MPI_INTEGER,&
status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, periodBsub, ndim-1, MPI_LOGICAL,&
status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, geometry_name(1:name_len), name_len,&
MPI_CHARACTER, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, stagger_grid, 1, MPI_LOGICAL, status,&
ierrmpi)
do iw = 1, nw
dname = trim(adjustl(cons_wnames(iw)))
call MPI_FILE_WRITE(slice_fh, dname, name_len, MPI_CHARACTER,&
status, ierrmpi)
end do
call MPI_FILE_WRITE(slice_fh, physics_type, name_len, MPI_CHARACTER,&
status, ierrmpi)
call phys_write_info(slice_fh)
if(pass_wall_time.or.save_now) then
call MPI_FILE_WRITE(slice_fh, snapshotnext, 1, MPI_INTEGER, status,&
ierrmpi)
else
call MPI_FILE_WRITE(slice_fh, snapshotnext+1, 1, MPI_INTEGER,&
status, ierrmpi)
end if
call MPI_FILE_WRITE(slice_fh, slicenext, 1, MPI_INTEGER, status,&
ierrmpi)
call MPI_FILE_WRITE(slice_fh, collapsenext, 1, MPI_INTEGER, status,&
ierrmpi)
call MPI_FILE_GET_POSITION(slice_fh, tree_offset, ierrmpi)
call MPI_FILE_WRITE(slice_fh, isleaf_sub, nsubleafs+nsubparents,&
MPI_LOGICAL, status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, block_lvl_sub, nsubleafs, MPI_INTEGER,&
status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, block_ig_sub, size(block_ig_sub),&
MPI_INTEGER, status, ierrmpi)
call MPI_FILE_GET_POSITION(slice_fh, block_offset, ierrmpi)
block_offset_sub(:) = block_offset_sub(:) + block_offset + &
int(nsubleafs*2*size_int, kind=MPI_OFFSET_KIND)
call MPI_FILE_WRITE(slice_fh, block_offset_sub, nsubleafs, MPI_OFFSET,&
status, ierrmpi)
call MPI_FILE_GET_POSITION(slice_fh, block_offset, ierrmpi)
if (block_offset - tree_offset /= &
(nsubleafs + nsubparents) * size_logical + &
nsubleafs * ((ndim-1+1) * size_int + 2 * size_int)) then
print *, "Warning: MPI_OFFSET type /= 8 bytes"
print *, "This *could* cause problems when reading .dat files"
end if
call MPI_FILE_SEEK(slice_fh, int(size_int,kind=MPI_OFFSET_KIND),&
MPI_SEEK_SET, ierrmpi)
call MPI_FILE_WRITE(slice_fh, int(tree_offset), 1, MPI_INTEGER,&
status, ierrmpi)
call MPI_FILE_WRITE(slice_fh, int(block_offset), 1, MPI_INTEGER,&
status, ierrmpi)
call MPI_FILE_CLOSE(slice_fh,ierrmpi)
end if
call MPI_BCAST(nsubleafs, 1, MPI_INTEGER, 0, icomm, ierrmpi)
call MPI_BCAST(nsubparents, 1, MPI_INTEGER, 0, icomm, ierrmpi)
if (mype /= 0) then
allocate(block_offset_sub(nsubleafs))
end if
call MPI_BCAST(block_offset_sub, nsubleafs, MPI_OFFSET, 0, icomm,&
ierrmpi)
amode=ior(MPI_MODE_CREATE,MPI_MODE_WRONLY)
call MPI_FILE_OPEN(icomm,filename,amode,MPI_INFO_NULL,slice_fh,ierrmpi)
iorequest=MPI_REQUEST_NULL
iwrite=0
n_ghost_sub = 0
do jgrid=1,Njgrid
iwrite=iwrite+1
offset=block_offset_sub(int(Morton_sub_start(mype)+jgrid-1))
call MPI_FILE_IWRITE_AT(slice_fh,offset,n_ghost_sub,2*(ndim-1),&
MPI_INTEGER, iorequest(iwrite),ierrmpi)
iwrite=iwrite+1
offset1=offset+int(2*(ndim-1)*size_int,kind=MPI_OFFSET_KIND)
call MPI_FILE_IWRITE_AT(slice_fh,offset1,ps_sub(jgrid)%w,1,&
type_subblock_io, iorequest(iwrite),ierrmpi)
end do
if (iwrite>0) call MPI_WAITALL(iwrite,iorequest,iostatus,ierrmpi)
call MPI_BARRIER(icomm, ierrmpi)
call MPI_FILE_CLOSE(slice_fh,ierrmpi)
deallocate(block_offset_sub)
if (mype==0) then
deallocate(block_ig_sub)
deallocate(block_lvl_sub)
deallocate(isleaf_sub)
end if
end subroutine put_slice_dat
subroutine put_slice_zerod
use mod_calculate_xw
integer:: iw
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
logical, save :: opened=.false.
character(len=1024) ::line, data
character(len=1024) :: filename, xlabel
character(len=79) :: xxlabel
integer :: amode, iwrite, status(MPI_STATUS_SIZE)
end subroutine put_slice_zerod
end subroutine put_slice