subroutine oneblock(qunit)
! this is for turning an AMR run into a single block
! the data will be all on selected level level_io
! this version should work for any dimension
! only writes w_write selected 1:nw variables, also nwauxio
! may use saveprim to switch to primitives
! this version can not work on multiple CPUs
! does not renormalize variables
! header info differs from onegrid below
! ASCII or binary output
use mod_forest, only: Morton_start, Morton_stop, sfc_to_igrid,&
igrid_to_node
use mod_global_parameters
use mod_usr_methods, only: usr_aux_output
use mod_physics
use mod_calculate_xw
integer, intent(in) :: qunit
integer :: Morton_no,igrid,ix1,ix2,ix3,ig1,ig2,ig3,level
integer, pointer :: ig_to_igrid(:,:,:,:)
logical :: fileopen,writeblk(max_blocks)
character(len=80) :: filename
integer :: filenr,ncells1,ncells2,ncells3,ncellx1,ncellx2,&
ncellx3,jg1,jg2,jg3,jig1,jig2,jig3
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
double precision :: wval1,xval1
double precision, dimension(1:1,1:1,1:1,1:nw+nwauxio) :: wval
double precision, dimension(1:1,1:1,1:1,1:ndim) :: xval
double precision:: normconv(0:nw+nwauxio)
integer :: iw,iiw,writenw,iwrite(1:nw+nwauxio),iigrid,idim
logical :: patchw(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3)
if(level_io<1)then
call mpistop('please specify level_io>0 for usage with oneblock')
end if
if(autoconvert)then
call mpistop('Set autoconvert=F and convert oneblock data manually')
end if
if(npe>1)then
if(mype==0) PRINT *,'ONEBLOCK as yet to be parallelized'
call mpistop('npe>1, oneblock')
end if
! only variables selected by w_write will be written out
normconv(0:nw+nwauxio)=one
normconv(0) = length_convert_factor
normconv(1:nw) = w_convert_factor
writenw=count(w_write(1:nw))+nwauxio
iiw=0
do iw =1,nw
if (.not.w_write(iw))cycle
iiw=iiw+1
iwrite(iiw)=iw
end do
if(nwauxio>0)then
do iw =nw+1,nw+nwauxio
iiw=iiw+1
iwrite(iiw)=iw
end do
end if
allocate(ig_to_igrid(ng1(level_io),ng2(level_io),ng3(level_io),0:npe-1))
ig_to_igrid=-1
writeblk=.false.
do Morton_no=Morton_start(mype),Morton_stop(mype)
igrid=sfc_to_igrid(Morton_no)
level=node(plevel_,igrid)
ig1=igrid_to_node(igrid,mype)%node%ig1
ig2=igrid_to_node(igrid,mype)%node%ig2
ig3=igrid_to_node(igrid,mype)%node%ig3;
ig_to_igrid(ig1,ig2,ig3,mype)=igrid
if((rnode(rpxmin1_,igrid)>=xprobmin1+&
(xprobmax1-xprobmin1)*writespshift(1,1).and.rnode(rpxmin2_,&
igrid)>=xprobmin2+(xprobmax2-xprobmin2)*writespshift(2,&
1).and.rnode(rpxmin3_,igrid)>=xprobmin3+&
(xprobmax3-xprobmin3)*writespshift(3,1)).and.(rnode(rpxmax1_,&
igrid)<=xprobmax1-(xprobmax1-xprobmin1)*writespshift(1,&
2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2-&
xprobmin2)*writespshift(2,2).and.rnode(rpxmax3_,&
igrid)<=xprobmax3-(xprobmax3-xprobmin3)*writespshift(3,2))) then
writeblk(igrid)=.true.
end if
end do
call getheadernames(wnamei,xandwnamei,outfilehead)
ncells1=0;ncells2=0;ncells3=0;
ncellx1=ixMhi1-ixMlo1+1
ncellx2=ixMhi2-ixMlo2+1
ncellx3=ixMhi3-ixMlo3+1
do ig1=1,ng1(level_io)
do ig2=1,ng2(level_io)
do ig3=1,ng3(level_io)
igrid=ig_to_igrid(ig1,ig2,ig3,mype)
if(writeblk(igrid)) go to 20
end do
end do
end do
20 continue
jg1=ig1;jg2=ig2;jg3=ig3;
jig1=jg1;jig2=jg2;jig3=jg3;
do ig1=1,ng1(level_io)
jig1=ig1
igrid=ig_to_igrid(jig1,jig2,jig3,mype)
if(writeblk(igrid)) ncells1=ncells1+ncellx1
end do
jig1=jg1;jig2=jg2;jig3=jg3;
do ig2=1,ng2(level_io)
jig2=ig2
igrid=ig_to_igrid(jig1,jig2,jig3,mype)
if(writeblk(igrid)) ncells2=ncells2+ncellx2
end do
jig1=jg1;jig2=jg2;jig3=jg3;
do ig3=1,ng3(level_io)
jig3=ig3
igrid=ig_to_igrid(jig1,jig2,jig3,mype)
if(writeblk(igrid)) ncells3=ncells3+ncellx3
end do
do iigrid=1,igridstail; igrid=igrids(iigrid)
if(.not.writeblk(igrid)) cycle
block=>ps(igrid)
if (nwauxio > 0) then
if (.not. associated(usr_aux_output)) then
call mpistop("usr_aux_output not defined")
else
call usr_aux_output(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
ixMlo1-1,ixMlo2-1,ixMlo3-1,ixMhi1+1,ixMhi2+1,ixMhi3+1,&
ps(igrid)%w,ps(igrid)%x,normconv)
end if
end if
end do
if (saveprim) then
do iigrid=1,igridstail; igrid=igrids(iigrid)
if (.not.writeblk(igrid)) cycle
block=>ps(igrid)
call phys_to_primitive(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
ixGlo1+1,ixGlo2+1,ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,ps(igrid)%w,&
ps(igrid)%x)
if(B0field) then
! add background magnetic field B0 to B
ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
iw_mag(:))=ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
iw_mag(:))+ps(igrid)%B0(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
:,0)
end if
end do
else
do iigrid=1,igridstail; igrid=igrids(iigrid)
if (.not.writeblk(igrid)) cycle
block=>ps(igrid)
if (B0field) then
! add background magnetic field B0 to B
if(phys_energy) ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
ixGlo3:ixGhi3,iw_e)=ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
ixGlo3:ixGhi3,iw_e)+0.5d0*sum(ps(igrid)%B0(ixGlo1:ixGhi1,&
ixGlo2:ixGhi2,ixGlo3:ixGhi3,:,0)**2,&
dim=ndim+1) + sum(ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
ixGlo3:ixGhi3,iw_mag(:))*ps(igrid)%B0(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
ixGlo3:ixGhi3,:,0),dim=ndim+1)
ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
iw_mag(:))=ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
iw_mag(:))+ps(igrid)%B0(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
:,0)
end if
end do
end if
Master_cpu_open : if (mype == 0) then
inquire(qunit,opened=fileopen)
if (.not.fileopen) then
! generate filename
filenr=snapshotini
if (autoconvert) filenr=snapshotnext
write(filename,'(a,i4.4,a)') TRIM(base_filename),filenr,".blk"
select case(convert_type)
case("oneblock")
open(qunit,file=filename,status='unknown')
write(qunit,*) TRIM(outfilehead)
write(qunit,*) ncells1,ncells2,ncells3
write(qunit,*) real(global_time*time_convert_factor)
case("oneblockB")
open(qunit,file=filename,form='unformatted',status='unknown')
write(qunit) outfilehead
write(qunit) ncells1,ncells2,ncells3
write(qunit) real(global_time*time_convert_factor)
end select
end if
end if Master_cpu_open
do ig3=1,ng3(level_io)
do ix3=ixMlo3,ixMhi3
do ig2=1,ng2(level_io)
do ix2=ixMlo2,ixMhi2
do ig1=1,ng1(level_io)
igrid=ig_to_igrid(ig1,ig2,ig3,mype)
if(.not.writeblk(igrid)) cycle
do ix1=ixMlo1,ixMhi1
Master_write : if(mype==0) then
select case(convert_type)
case("oneblock")
write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix1,ix2,ix3,&
1:ndim)*normconv(0),(ps(igrid)%w(ix1,ix2,ix3,&
iwrite(iw))*normconv(iwrite(iw)),iw=1,writenw)
case("oneblockB")
write(qunit) real(ps(igrid)%x(ix1,ix2,ix3,&
1:ndim)*normconv(0)),(real(ps(igrid)%w(ix1,ix2,ix3,&
iwrite(iw))*normconv(iwrite(iw))),iw=1,writenw)
end select
end if Master_write
end do
end do
end do
end do
end do
end do
close(qunit)
end subroutine oneblock