onegrid Subroutine

  subroutine onegrid(qunit)
    ! this is for turning an AMR run into a single grid
    ! this version should work for any dimension, can be in parallel
    ! in 1D, should behave much like oneblock, except for header info

    ! only writes all 1:nw variables, no nwauxio
    ! may use saveprim to switch to primitives
    ! this version can work on multiple CPUs
    ! does not renormalize variables
    ! ASCII output

    use mod_forest, only: Morton_start, Morton_stop, sfc_to_igrid
    use mod_global_parameters
    use mod_physics
    use mod_calculate_xw
    integer, intent(in) :: qunit

    integer             :: itag,Morton_no,igrid,ix1,ix2,ix3,iw
    logical             :: fileopen
    character(len=80)   :: filename
    integer             :: filenr

    character(len=name_len) :: wnamei(1:nw+nwauxio),&
       xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead

    !.. MPI variables ..
    integer           :: igrid_recv,ipe
    double precision  :: w_recv(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:nw),x_recv(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:ndim)
    integer, allocatable :: intstatus(:,:)

    if(nwauxio>0)then
      if(mype==0) PRINT *,'ONEGRID to be used without nwauxio'
      call mpistop('nwauxio>0, onegrid')
    end if

    if(saveprim)then
      if(mype==0.and.nwaux>0) PRINT *&
         ,'warning: ONEGRID used with saveprim, check auxiliaries'
    end if

    Master_cpu_open : if (mype == 0) then
      call getheadernames(wnamei,xandwnamei,outfilehead)
      write(outfilehead,'(a)') "#"//" "//TRIM(outfilehead)
      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"
        open(qunit,file=filename,status='unknown')
      end if
      write(qunit,"(a)")outfilehead
      write(qunit,"(i7)") ( (ixMhi1-ixMlo1+1)*(ixMhi2-ixMlo2+1)*(ixMhi3-ixMlo3+&
         1) )*(Morton_stop(npe-1)-Morton_start(0)+1)
    end if Master_cpu_open

    do Morton_no=Morton_start(mype),Morton_stop(mype)
      igrid=sfc_to_igrid(Morton_no)
      if(saveprim) call phys_to_primitive(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,&
         ixGhi3,ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,ps(igrid)%w,&
         ps(igrid)%x)
      if(mype/=0)then
        itag=Morton_no
        call mpi_send_wrapper(igrid,1,MPI_INTEGER, 0,itag,icomm,ierrmpi)
        call mpi_send_wrapper(ps(igrid)%x,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
        itag=igrid
        call mpi_send_wrapper(ps(igrid)%w,1,type_block_io, 0,itag,icomm,ierrmpi)
      else
       do ix3=ixMlo3,ixMhi3
       do ix2=ixMlo2,ixMhi2
       do ix1=ixMlo1,ixMhi1
          do iw=1,nw
            if( dabs(ps(igrid)%w(ix1,ix2,ix3,iw)) < 1.0d-32 ) ps(igrid)%w(ix1,&
               ix2,ix3,iw) = zero
          end do
          write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix1,ix2,ix3,1:ndim),&
             ps(igrid)%w(ix1,ix2,ix3,1:nw)
       end do
       end do
       end do
      end if
    end do

    if(mype==0.and.npe>1) allocate(intstatus(MPI_STATUS_SIZE,1))

    Manycpu : if (npe>1) then
      if (mype==0) then
       loop_cpu : do ipe =1, npe-1
        loop_Morton : do Morton_no=Morton_start(ipe),Morton_stop(ipe)
              itag=Morton_no
              call mpi_recv_wrapper(igrid_recv,1,MPI_INTEGER, ipe,itag,icomm,&
                 intstatus(:,1),ierrmpi)
              call mpi_recv_wrapper(x_recv,1,type_block_xcc_io, ipe,itag,icomm,&
                 intstatus(:,1),ierrmpi)
              itag=igrid_recv
              call mpi_recv_wrapper(w_recv,1,type_block_io, ipe,itag,icomm,intstatus(:,&
                 1),ierrmpi)
              do ix3=ixMlo3,ixMhi3
              do ix2=ixMlo2,ixMhi2
              do ix1=ixMlo1,ixMhi1
                 do iw=1,nw
                   if( dabs(ps(igrid)%w(ix1,ix2,ix3,&
                      iw)) < smalldouble ) ps(igrid)%w(ix1,ix2,ix3,iw) = zero
                 end do
                 write(qunit,fmt="(100(e14.6))") x_recv(ix1,ix2,ix3,1:ndim),&
                    w_recv(ix1,ix2,ix3,1:nw)
              end do
              end do
              end do
        end do loop_Morton
       end do loop_cpu
      end if
    end if Manycpu

    if (npe>1) then
      call MPI_BARRIER(icomm,ierrmpi)
      if(mype==0)deallocate(intstatus)
    end if

    if(mype==0) close(qunit)
  end subroutine onegrid