mod_collapse.t

Go to the documentation of this file.

!> Collapses D-dimensional output to D-1 view by line-of-sight integration
! for now: only for non-stretched cartesian cases, with LOS along coordinate
! writes out as either csv or vti file (collapse_type)
module mod_collapse
  use mod_global_parameters
  implicit none
 
contains
!=============================================================================
subroutine write_collapsed
use mod_global_parameters
! Writes a collapsed view of the data integrated over one grid-direction.  
! E.g. column density maps.  
! Uses flat interpolation throughout.
! by Oliver Porth
! 6.Nov 2013
integer :: idir
logical, save :: firstcollapse=.true.
!-----------------------------------------------------------------------------
 
if (firstcollapse) then
   firstcollapse=.false.
end if
 
do idir=1, ndim
   if (collapse(idir)) call put_collapse(idir)
end do
 
collapsenext=collapsenext+1
end subroutine write_collapsed
!=============================================================================
subroutine put_collapse(dir)
use mod_forest, only: morton_start, morton_stop, sfc_to_igrid
use mod_slice, only:alloc_subnode, dealloc_subnode
use mod_global_parameters
integer, intent(in)                               :: dir
! .. local ..
integer                                           :: jgrid, igrid, Morton_no
double precision,dimension(0:nw+nwauxio)          :: normconv 
!-----------------------------------------------------------------------------
if(.not.slab) call mpistop("collapse only for slab cartesian cases")
 
call allocate_collapsed(dir)
 
jgrid=0
do morton_no=morton_start(mype),morton_stop(mype)
   igrid=sfc_to_igrid(morton_no)
   jgrid=jgrid+1
   call alloc_subnode(jgrid,dir,nwauxio)
   call collapse_subnode(igrid,jgrid,dir,normconv)
   call integrate_subnode(igrid,jgrid,dir)
end do
 
! Reduce to head-node:
if (mype==0) then
   call mpi_reduce(mpi_in_place,collapseddata,size(collapseddata),mpi_double_precision,mpi_sum,0,icomm,ierrmpi)
else
   call mpi_reduce(collapseddata,collapseddata,size(collapseddata),mpi_double_precision,mpi_sum,0,icomm,ierrmpi)
end if
call mpi_barrier(icomm, ierrmpi)
 
{^nooned
select case(collapse_type)
case('vti')
   call output_collapsed_vti(dir,normconv)
case('csv')
   call output_collapsed_csv(dir,normconv)
case('default')
   call mpistop("Unknown filetype for collapsed views")
end select
   }{^ifoned
call mpistop("sorry, 1D collapsed output routine not yet implemented (should be easy)...")
}
 
! If we need the subnodes later, remove deallocation here:
do jgrid=1,morton_stop(mype)-morton_start(mype)+1
   call dealloc_subnode(jgrid)
end do
deallocate(collapseddata)
 
end subroutine put_collapse
!=============================================================================
subroutine output_collapsed_csv(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
{^nooned
integer                                           :: ix^DM
double precision                                  :: dxdim^DM, xdim^LIM^DM
}
!-----------------------------------------------------------------------------
 
if (mype==0) then
 
! Get coordinates:
{^ifthreed
select case(dir)
case (1)
   dxdim1 = dx(2,collapselevel)
   dxdim2 = dx(3,collapselevel)
   xdim^lim1=xprob^lim2; 
   xdim^lim2=xprob^lim3; 
case (2)
   dxdim1 = dx(1,collapselevel)
   dxdim2 = dx(3,collapselevel)
   xdim^lim1=xprob^lim1; 
   xdim^lim2=xprob^lim3; 
case (3)
   dxdim1 = dx(1,collapselevel)
   dxdim2 = dx(2,collapselevel)
   xdim^lim1=xprob^lim1; 
   xdim^lim2=xprob^lim2; 
case default
   dxdim1 = 1
   dxdim2 = 1
   xdim^lim1=xprob^lim2; 
   xdim^lim2=xprob^lim3; 
   call mpistop("slice direction not clear in output_collapsed_csv")
end select
}
{^iftwod
select case(dir)
case (1)
   dxdim1 = dx(2,collapselevel)
   xdim^lim1=xprob^lim2; 
case (2)
   dxdim1 = dx(1,collapselevel)
   xdim^lim1=xprob^lim1; 
case default
   dxdim1 = 1
   xdim^lim1=xprob^lim2; 
   call mpistop("slice direction not clear in output_collapsed_csv")
end select
}
 
 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,'.csv'
      open(unitcollapse,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
      if (iw .eq. dir) cycle
      line = trim(line)//trim(xandwnamei(iw))//', '
   end do
   line = trim(line)//trim(xandwnamei(ndim+nw+nwauxio))
   write(unitcollapse,'(a)')trim(line)
   myshape = shape(collapseddata)
{^nooned
{^dm& do ix^dmb = 1,myshape(^dmb)\}
   ! following assumes uniform cartesian grid
   write(unitcollapse,'(200(1pe20.12))') &
        {^dm& dxdim^dm*dble(ix^dm)+xdimmin^dm}, &
        (collapseddata(ix^dm,iw)*normconv(iw),iw=1,nw+nwauxio)
{^dm& enddo\}
}
close(unitcollapse)
 
end if
end subroutine output_collapsed_csv
!=============================================================================
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
{^nooned
integer                                           :: ix^DM
double precision                                  :: dxdim^DM, xdim^LIM^DM
}
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:
{^ifthreed
select case(dir)
case (1)
   dxdim1 = dx(2,1)*2.0d0**(1-collapselevel)
   dxdim2 = dx(3,1)*2.0d0**(1-collapselevel)
   xdim^lim1=xprob^lim2; 
   xdim^lim2=xprob^lim3; 
case (2)
   dxdim1 = dx(1,1)*2.0d0**(1-collapselevel)
   dxdim2 = dx(3,1)*2.0d0**(1-collapselevel)
   xdim^lim1=xprob^lim1; 
   xdim^lim2=xprob^lim3; 
case (3)
   dxdim1 = dx(1,1)*2.0d0**(1-collapselevel)
   dxdim2 = dx(2,1)*2.0d0**(1-collapselevel)
   xdim^lim1=xprob^lim1; 
   xdim^lim2=xprob^lim2; 
case default
   dxdim1 = 1
   dxdim2 = 1
   xdim^lim1=xprob^lim2; 
   xdim^lim2=xprob^lim3; 
   call mpistop("slice direction not clear in output_collapsed_vti")
end select
}
{^iftwod
select case(dir)
case (1)
   dxdim1 = dx(2,1)*2.0d0**(1-collapselevel)
   xdim^lim1=xprob^lim2; 
case (2)
   dxdim1 = dx(1,1)*2.0d0**(1-collapselevel)
   xdim^lim1=xprob^lim1; 
case default
   dxdim1 = 1
   xdim^lim1=xprob^lim2; 
   call mpistop("slice direction not clear in output_collapsed_vti")
end select
}
 
origin=0
spacing=zero
wholeextent=0
myshape = shape(collapseddata)
{^ifoned
length = size_single
}
{^nooned
length = ^dm&myshape(^dm)*
length = length*size_single
{^dm&wholeextent(^dm*2)=myshape(^dm); }
{^dm&spacing(^dm) = dxdim^dm; }
{^dm&origin(^d) = xdimmin^dm; }
}
 
 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
  {^ifoned
  write(unitcollapse) real(collapseddata(iw)*normconv(iw))
  }
  {^iftwod
  write(unitcollapse) (real(collapseddata(ix1,iw)*normconv(iw)),ix1=1,myshape(1))
  }
  {^ifthreed
  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
!=============================================================================
subroutine allocate_collapsed(dir)
use mod_global_parameters
integer, intent(in)                               :: dir
integer                                           :: dim^D
integer                                           :: nx^D
!-----------------------------------------------------------------------------
! allocate array for the collapsed data:
! number of cells per grid.
nx^d=ixmhi^d-ixmlo^d+1;
{dim^d=ng^d(1)*2**(collapselevel-1)*nx^d\}
{^ifthreed
select case(dir)
case (1)
   allocate(collapseddata(&
        dim2,dim3,nw+nwauxio))
case (2)
   allocate(collapseddata(&
        dim1,dim3,nw+nwauxio))
case (3)
   allocate(collapseddata(&
        dim1,dim2,nw+nwauxio))
case default
   call mpistop("slice direction not clear in allocate_collapsed")
end select
}
{^iftwod
select case(dir)
case (1)
   allocate(collapseddata(&
        dim2,nw+nwauxio))
case (2)
   allocate(collapseddata(&
        dim1,nw+nwauxio))
case default
   call mpistop("slice direction not clear in allocate_collapsed")
end select
}
{^ifoned
   allocate(collapseddata(nw+nwauxio))
}
collapseddata = zero
 
end subroutine allocate_collapsed
!=============================================================================
subroutine integrate_subnode(igrid,jgrid,dir)
use mod_forest, only: tree_node_ptr, igrid_to_node
use mod_global_parameters
integer, intent(in)                        :: igrid, jgrid, dir
! .. local ..
type(tree_node_ptr)                        :: tree
integer                                    :: nx^D
integer                                    :: ig^D, level, ig^Dtargetmin, ig^Dtargetmax
integer                                    :: ix^Dtarget^LIM, idim^Dtarget^LIM
integer                                    :: ixMdim^LLIM^D, ix^Dorig, ix^D
{^nooned
integer                                    ::  igdim^DM
}
!-----------------------------------------------------------------------------
tree%node => igrid_to_node(igrid, mype)%node
{^d& ig^d = tree%node%ig^d; }
level = tree%node%level
! number of cells per grid.
nx^d=ixmhi^d-ixmlo^d+1;
 
! assumes uniform cartesian grid with factor 2 refinement per dimension
{^d&
if (level > collapselevel) then
   ig^dtargetmin = int(dble(ig^d-1)/2.0d0**(level-collapselevel))+1
   ig^dtargetmax = ig^dtargetmin
else if (level < collapselevel) then
   ig^dtargetmin = int(2.0d0**(collapselevel-level))*(ig^d-1)+1
   ig^dtargetmax = int(2.0d0**(collapselevel-level))*ig^d
else
   ig^dtargetmin = ig^d
   ig^dtargetmax = ig^d
end if
ix^dtargetmin = nx^d*(ig^dtargetmin-1)+1
ix^dtargetmax = nx^d*ig^dtargetmax
\}
 
{^ifthreed
select case(dir)
case (1)
   igdim1=(ig2-1)*nx2
   igdim2=(ig3-1)*nx3
   idim1target^lim=ix2target^lim; 
   idim2target^lim=ix3target^lim;
   ixmdim^llim1=ixm^llim2;
   ixmdim^llim2=ixm^llim3;
case (2)
   igdim1=(ig1-1)*nx1
   igdim2=(ig3-1)*nx3
   idim1target^lim=ix1target^lim; 
   idim2target^lim=ix3target^lim;
   ixmdim^llim1=ixm^llim1;
   ixmdim^llim2=ixm^llim3;
case (3)
   igdim1=(ig1-1)*nx1
   igdim2=(ig2-1)*nx2
   idim1target^lim=ix1target^lim; 
   idim2target^lim=ix2target^lim;
   ixmdim^llim1=ixm^llim1;
   ixmdim^llim2=ixm^llim2;
case default
   call mpistop("slice direction not clear in integrate_subnode")
end select
 
if (level >= collapselevel) then
   do ix1orig = ixmdimlo1,ixmdimhi1
      do ix2orig = ixmdimlo2,ixmdimhi2
{^dm& ix^dm = int(dble(ix^dmorig-nghostcells+igdim^dm-1)*2.0d0**(collapselevel-level))+1\}
         collapseddata(ix1,ix2,1:nw+nwauxio) = collapseddata(ix1,ix2,1:nw+nwauxio) &
              + ps_sub(jgrid)%w(ix1orig,ix2orig,1:nw+nwauxio) / 2.0d0**(2*(level-collapselevel))
      end do
   end do
else
{^dm&
   do ix^dm = idim^dmtargetmin,idim^dmtargetmax\}
 {^dm& ix^dmorig = int(dble(ix^dm-idim^dmtargetmin)/2.0d0**(collapselevel-level))+1+nghostcells \}
 collapseddata(ix1,ix2,1:nw+nwauxio) = collapseddata(ix1,ix2,1:nw+nwauxio) + ps_sub(jgrid)%w(ix1orig,ix2orig,1:nw+nwauxio)
 {^dm& enddo\}
end if
}
{^iftwod
select case(dir)
case (1)
   igdim1=(ig2-1)*nx2
   idim1target^lim=ix2target^lim; 
   ixmdim^llim1=ixm^llim2;
case (2)
   igdim1=(ig1-1)*nx1
   idim1target^lim=ix1target^lim; 
   ixmdim^llim1=ixm^llim1;
case default
   call mpistop("slice direction not clear in integrate_subnode")
end select
 
if (level >= collapselevel) then
   do ix1orig = ixmdimlo1,ixmdimhi1
{^dm& ix^dm = int(dble(ix^dmorig-nghostcells+igdim^dm-1)*2.0d0**(collapselevel-level))+1\}
         collapseddata(ix1,1:nw+nwauxio) = collapseddata(ix1,1:nw+nwauxio) &
              + ps_sub(jgrid)%w(ix1orig,1:nw+nwauxio) / 2.0d0**(level-collapselevel)
   end do
else
{^dm&
   do ix^dm = idim^dmtargetmin,idim^dmtargetmax\}
 {^dm& ix^dmorig = int(dble(ix^dm-idim^dmtargetmin)/2.0d0**(collapselevel-level))+1+nghostcells \}
 collapseddata(ix1,1:nw+nwauxio) = collapseddata(ix1,1:nw+nwauxio) + ps_sub(jgrid)%w(ix1orig,1:nw+nwauxio)
 {^dm& enddo\}
end if
}
{^ifoned
collapseddata(1:nw+nwauxio) = collapseddata(1:nw+nwauxio) + ps_sub(jgrid)%w(1:nw+nwauxio)
}
 
end subroutine integrate_subnode
!=============================================================================
subroutine collapse_subnode(igrid,jgrid,dir,normconv)
  use mod_usr_methods, only: usr_aux_output
  use mod_global_parameters
  use mod_physics, only: phys_to_primitive
  use mod_calculate_xw
 
integer, intent(in) :: igrid, jgrid, dir
double precision,dimension(0:nw+nwauxio),intent(out)       :: normconv 
! .. local ..
integer                :: ix, iw
double precision       :: dx^D
double precision, dimension(ixMlo^D-1:ixMhi^D,ndim)       :: xC_TMP, xC
double precision, dimension(ixMlo^D:ixMhi^D,ndim)         :: xCC_TMP, xCC
double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)   :: wCC_TMP
integer                :: ixC^L, ixCC^L
!-----------------------------------------------------------------------------
ps_sub(jgrid)%w=zero
dx^d=rnode(rpdx^d_,igrid);
 
call calc_x(igrid,xc,xcc)
call calc_grid(unitslice,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
         ixc^l,ixcc^l,.true.)
 
{^ifthreed
select case (dir)
case (1)
  if(slab_uniform) then
    do ix=ixmlo1,ixmhi1
      ps_sub(jgrid)%w(ixmlo2:ixmhi2,ixmlo3:ixmhi3,1:nw+nwauxio) = &
           ps_sub(jgrid)%w(ixmlo2:ixmhi2,ixmlo3:ixmhi3,1:nw+nwauxio) &
           + wcc_tmp(ix,ixmlo2:ixmhi2,ixmlo3:ixmhi3,1:nw+nwauxio) * dx1
    end do
  else
    do iw=1,nw+nwauxio
      do ix=ixmlo1,ixmhi1
        ps_sub(jgrid)%w(ixmlo2:ixmhi2,ixmlo3:ixmhi3,iw) = &
             ps_sub(jgrid)%w(ixmlo2:ixmhi2,ixmlo3:ixmhi3,iw) &
             + wcc_tmp(ix,ixmlo2:ixmhi2,ixmlo3:ixmhi3,iw) * &
             ps(igrid)%dx(ix,ixmlo2:ixmhi2,ixmlo3:ixmhi3,1)
      end do
    end do
  end if
  ps_sub(jgrid)%x(ixmlo2:ixmhi2,ixmlo3:ixmhi3,1:ndim) = &
       ps(igrid)%x(ixmlo1,ixmlo2:ixmhi2,ixmlo3:ixmhi3,1:ndim)
case (2)
  if(slab_uniform) then
    do ix=ixmlo2,ixmhi2
       ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo3:ixmhi3,1:nw+nwauxio) = &
            ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo3:ixmhi3,1:nw+nwauxio) &
            + wcc_tmp(ixmlo1:ixmhi1,ix,ixmlo3:ixmhi3,1:nw+nwauxio) * dx2
    end do
  else
    do iw=1,nw+nwauxio
      do ix=ixmlo2,ixmhi2
         ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo3:ixmhi3,iw) = &
              ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo3:ixmhi3,iw) &
              + wcc_tmp(ixmlo1:ixmhi1,ix,ixmlo3:ixmhi3,iw) * &
             ps(igrid)%dx(ixmlo1:ixmhi1,ix,ixmlo3:ixmhi3,2)
      end do
    end do
  endif
  ps_sub(jgrid)%x(ixmlo1:ixmhi1,ixmlo3:ixmhi3,1:ndim) = &
       ps(igrid)%x(ixmlo1:ixmhi1,ixmlo2,ixmlo3:ixmhi3,1:ndim) 
case (3)
  if(slab_uniform) then
    do ix=ixmlo3,ixmhi3
       ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo2:ixmhi2,1:nw+nwauxio) = &
            ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo2:ixmhi2,1:nw+nwauxio) &
            + wcc_tmp(ixmlo1:ixmhi1,ixmlo2:ixmhi2,ix,1:nw+nwauxio) * dx3
    end do
  else
    do iw=1,nw+nwauxio
      do ix=ixmlo3,ixmhi3
         ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo2:ixmhi2,iw) = &
              ps_sub(jgrid)%w(ixmlo1:ixmhi1,ixmlo2:ixmhi2,iw) &
              + wcc_tmp(ixmlo1:ixmhi1,ixmlo2:ixmhi2,ix,iw) * &
             ps(igrid)%dx(ixmlo1:ixmhi1,ixmlo2:ixmhi2,ix,3)
      end do
    end do
  endif
  ps_sub(jgrid)%x(ixmlo1:ixmhi1,ixmlo2:ixmhi2,1:ndim) = &
       ps(igrid)%x(ixmlo1:ixmhi1,ixmlo2:ixmhi2,ixmlo3,1:ndim) 
case default
   print*, 'subnode, dir: ', dir
   call mpistop("slice direction not clear in collapse_subnode")
end select
}
{^iftwod
select case (dir)
case (1)
  if(slab_uniform) then
    do ix=ixmlo1,ixmhi1
      ps_sub(jgrid)%w(ixmlo2:ixmhi2,1:nw+nwauxio) = &
           ps_sub(jgrid)%w(ixmlo2:ixmhi2,1:nw+nwauxio) &
           + wcc_tmp(ix,ixmlo2:ixmhi2,1:nw+nwauxio) * dx1
    end do
  else
    do iw=1,nw+nwauxio
      do ix=ixmlo1,ixmhi1
        ps_sub(jgrid)%w(ixmlo2:ixmhi2,iw) = &
             ps_sub(jgrid)%w(ixmlo2:ixmhi2,iw) &
             + wcc_tmp(ix,ixmlo2:ixmhi2,iw) * &
              ps(igrid)%dx(ix,ixmlo2:ixmhi2,1)
      end do
    end do
  end if
  ps_sub(jgrid)%x(ixmlo2:ixmhi2,1:ndim) = &
       ps(igrid)%x(ixmlo1,ixmlo2:ixmhi2,1:ndim)
case (2)
  if(slab_uniform) then
    do ix=ixmlo2,ixmhi2
       ps_sub(jgrid)%w(ixmlo1:ixmhi1,1:nw+nwauxio) = &
            ps_sub(jgrid)%w(ixmlo1:ixmhi1,1:nw+nwauxio) &
            + wcc_tmp(ixmlo1:ixmhi1,ix,1:nw+nwauxio) * dx2
    end do
  else
    do iw=1,nw+nwauxio
      do ix=ixmlo2,ixmhi2
         ps_sub(jgrid)%w(ixmlo1:ixmhi1,iw) = &
              ps_sub(jgrid)%w(ixmlo1:ixmhi1,iw) &
              + wcc_tmp(ixmlo1:ixmhi1,ix,iw) * &
              ps(igrid)%dx(ixmlo1:ixmhi1,ix,2)
      end do
    end do
  end if
  ps_sub(jgrid)%x(ixmlo1:ixmhi1,1:ndim) = &
       ps(igrid)%x(ixmlo1:ixmhi1,ixmlo2,1:ndim) 
case default
   call mpistop("slice direction not clear in collapse_subnode")
end select
}
{^ifoned   
if(slab_uniform) then
  do ix=ixmlo1,ixmhi1
     ps_sub(jgrid)%w(1:nw+nwauxio) = ps_sub(jgrid)%w(1:nw+nwauxio) + wcc_tmp(ix,1:nw+nwauxio) * dx1
  end do
else
  do iw=1,nw+nwauxio
    do ix=ixmlo1,ixmhi1
      ps_sub(jgrid)%w(iw) = ps_sub(jgrid)%w(iw) + wcc_tmp(ix,iw) * ps(igrid)%dx(ix,1)
    end do
  end do
end if
ps_sub(jgrid)%x(1:ndim) = ps(igrid)%x(ixmlo1,1:ndim)
}
 
end subroutine collapse_subnode
!=============================================================================
end module mod_collapse