Module class_DataDecomposition    !============================================
 use class_TensorProductGrid
 implicit none
 public  :: NewDataDecomposition
 public  :: DeleteDataDecomposition
 public  :: DescribeDecomposition
 
 type DataDecomposition 
  integer                            :: level         !!multigrid level 0=fine grid
  integer                            :: nx            !!global data size, dim 1
  integer                            :: ny            !!global data size, dim 2
  integer                            :: nz            !!global data size, dim 3
  integer                            :: numprocs            !!number of processors
  integer,pointer,dimension(:,:)     :: SlabDims,ColDims    !!size of local storage arrays
  integer,pointer,dimension(:,:)     :: SlabVals,ColVals    !!num of data vals held 
  integer,pointer,dimension(:)       :: my_global_kvals     !!global indexing of local data
  integer,pointer,dimension(:)       :: my_global_ivals     !!global indexing of local data
  integer                            :: locnz
  integer                            :: locnx
  character (len=80 )                :: DerivTypes(6)
  character (len=80 )                :: label
 end type DataDecomposition
 
 interface assignment (=)
  module procedure EqualDataDecompositions
 end interface
 
Contains         

 subroutine DeleteDataDecomposition(DataMap)
 type( DataDecomposition ),intent(inout)  :: DataMap
  deallocate( DataMap%SlabDims )
  deallocate( DataMap%ColDims )
  deallocate( DataMap%SlabVals )
  deallocate( DataMap%ColVals )
  deallocate( DataMap%my_global_kvals )
  deallocate( DataMap%my_global_ivals )
 end subroutine DeleteDataDecomposition
 
 
 subroutine EqualDataDecompositions(new,old)
  implicit none
  type( DataDecomposition ),intent(inout):: new
  type( DataDecomposition ),intent(in)   :: old
  integer                                :: numprocs,rcode(6)=0
  
  numprocs=old%numprocs
  
   if( associated(new%SlabDims) ) call DeleteDataDecomposition(new)
   
   new%level = old%level
   new%nx = old%nx
   new%ny = old%ny
   new%nz = old%nz
   new%numprocs = old%numprocs
   new%DerivTypes = old%DerivTypes
   new%label = old%label
   new%locnz = old%locnz
   new%locnx = old%locnx
   
   allocate( new%SlabDims(3,-1:numprocs-1), stat=rcode(1) )
   allocate( new%SlabVals(3,-1:numprocs-1), stat=rcode(2) )
   allocate( new%ColDims(3,-1:numprocs-1), stat=rcode(3) )
   allocate( new%ColVals(3,-1:numprocs-1), stat=rcode(4) )
   allocate( new%my_global_kvals(new%locnz), stat=rcode(5) ) 
   allocate( new%my_global_ivals(new%locnx), stat=rcode(6) ) 
   if( any( rcode /= 0 ) ) Stop 'Allocation failed: EqualDataDecompositions'   
     
   new%SlabDims = old%SlabDims
   new%ColDims = old%ColDims
   new%SlabVals = old%SlabVals
   new%ColVals = old%ColVals  
   new%my_global_kvals = old%my_global_kvals
   new%my_global_ivals = old%my_global_ivals
 end subroutine EqualDataDecompositions
 
 subroutine NewDataDecomposition(label,grid,level,DerivTypes,DataMap)
  use class_TensorProductGrid
  implicit none
   
  type( DataDecomposition ),intent(inout)      :: DataMap
  type( TensorProductGrid ),intent(in)         :: grid
  character (len=80 )                          :: DerivTypes(6)
  character (len=80 ),intent(in)               :: label
  integer                                      :: level
  integer                                      :: numnodes
  integer                                      :: i,k0,rem
  integer                                      :: myid,ierr,rcode(6)=0
  include 'mpif.h'
  
  call mpi_comm_rank (MPI_COMM_WORLD,myid,ierr)     
  call mpi_comm_size (MPI_COMM_WORLD,numnodes,ierr) 
  
  rem = mod(DataMap%nx,numnodes)
  if( rem /= 0 .and. rem /= 1 ) then
   if(myid==0) &
   write(0,*) ' nx should be evenly divisible by np or have remainder 1'
   stop
  endif
  
  rem = mod(DataMap%nz,numnodes)
  if( rem /= 0 .and. rem /= 1 ) then
   if(myid==0) &
   write(0,*) ' nz should be evenly divisible by np or have remainder 1'
   stop
  endif
  
  DataMap%level=level
  DataMap%label=label
  DataMap%DerivTypes(1:3) = DerivTypes(1:3)
  DataMap%nx=grid%DataDims(1)
  DataMap%ny=grid%DataDims(2)
  DataMap%nz=grid%DataDims(3)
  DataMap%numprocs=numnodes
  
  do i=1,3
   if( trim( DerivTypes(i) ) == 'cos' ) then
    DataMap%DerivTypes(i+3)='sin'
   elseif( trim( DerivTypes(i) ) == 'sin' ) then
    DataMap%DerivTypes(i+3)='cos'
   else
    DataMap%DerivTypes(i+3)=DerivTypes(i)
   endif
  enddo
 
  if( .not. associated(DataMap%SlabDims) ) then
   allocate( DataMap%SlabDims(3,-1:numnodes-1), stat=rcode(1) )
   allocate( DataMap%SlabVals(3,-1:numnodes-1), stat=rcode(2) )
   allocate( DataMap%ColDims(3,-1:numnodes-1), stat=rcode(3) )
   allocate( DataMap%ColVals(3,-1:numnodes-1), stat=rcode(4) )
   if( any( rcode /= 0 ) ) Stop 'Allocation failed: NewDataDecomposition A'
  endif
  
  i=1    
   if( DerivTypes(i) == 'fourier' ) then
    DataMap%SlabVals(i,:) = DataMap%nx
    DataMap%SlabDims(i,:) = DataMap%nx + 0        !! DON"T accomodate conjugate symmetry
    DataMap%ColDims(i,:)  = (DataMap%nx)/numnodes !! size counted in real words, but the
    DataMap%ColVals(i,:)  = (DataMap%nx)/numnodes !! storage is for local part of kx wavenumber space
   else
    select case ( mod(DataMap%nx,2) )  
     case(0)              ! nx is even
      DataMap%SlabVals(i,:) = DataMap%nx
      DataMap%SlabDims(i,:) = DataMap%nx
      DataMap%ColDims(i,:)  = DataMap%nx/numnodes
      DataMap%ColVals(i,:)  = DataMap%nx/numnodes
 
     case(1)              ! nx is odd
      DataMap%SlabVals(i,:) = DataMap%nx
      DataMap%SlabDims(i,:) = DataMap%nx
      DataMap%ColDims(i,:)  = (DataMap%nx-1)/numnodes
      DataMap%ColVals(i,:)  = (DataMap%nx-1)/numnodes
      DataMap%ColDims(i,numnodes-1) = (DataMap%nx-1)/numnodes + 1
      DataMap%ColVals(i,numnodes-1) = (DataMap%nx-1)/numnodes + 1
     end select
   endif
   
   i=2    
   DataMap%SlabVals(i,:) = DataMap%ny
   DataMap%SlabDims(i,:) = DataMap%ny
   DataMap%ColDims(i,:)  = DataMap%ny
   DataMap%ColVals(i,:)  = DataMap%ny
       
   i=3
    select case ( mod(DataMap%nz,2) )  
     case(0)              ! nz is even
      DataMap%SlabVals(i,:) = DataMap%nz/numnodes
      DataMap%SlabDims(i,:) = DataMap%nz/numnodes
      DataMap%ColDims(i,:) = DataMap%nz
      DataMap%ColVals(i,:) = DataMap%nz
     case(1)              ! nz is odd
      DataMap%SlabVals(i,:) = (DataMap%nz-1)/numnodes
      DataMap%SlabDims(i,:) = (DataMap%nz-1)/numnodes
      DataMap%ColDims(i,:)  = DataMap%nz
      DataMap%ColVals(i,:)  = DataMap%nz
      DataMap%SlabDims(i,numnodes-1) = (DataMap%nz-1)/numnodes + 1
      DataMap%SlabVals(i,numnodes-1) = (DataMap%nz-1)/numnodes + 1
     end select
     
   !! determine the global k indices for processor myid
   !! first, I already know how many planes are allocated to myid
     DataMap%locnz = DataMap%SlabDims(3,myid)
     
   !! next, add up the k locations below myid
     k0=0
     do i=0,myid-1  
      k0=k0+DataMap%SlabDims(3,i)  !! num of k planes assigned to processor i
     enddo     
     
     allocate(DataMap%my_global_kvals( DataMap%locnz ), stat=rcode(1) )
     if( rcode(1) /= 0  ) Stop 'Allocation failed: NewDataDecomposition B'
   !!fill the array w/ the k values
     do i=1,DataMap%locnz
       DataMap%my_global_kvals(i) = k0+i
     enddo
    
   !! determine the global i indices for processor myid (Column orientation)
   !! first, I already know how many planes are allocated to myid
     DataMap%locnx = DataMap%ColDims(1,myid)
     
   !! next, add up the i locations left of myid
     k0=0
     do i=0,myid-1  
      k0=k0+DataMap%ColDims(1,i)  !! num of X planes assigned to processor i
     enddo     
     
     allocate(DataMap%my_global_ivals( DataMap%locnx ), stat=rcode(1) )
     if( rcode(1) /= 0  ) Stop 'Allocation failed: NewDataDecomposition C'
   !!fill the array w/ the i values
     do i=1,DataMap%locnx
       DataMap%my_global_ivals(i) = k0+i
     enddo

    if(myid==0) write(0,*) '......creating MPI derived data types'
    call create_MPI_data_types(DataMap)
    if(myid==0) write(0,*) '....... data types done'
    
  end subroutine NewDataDecomposition
 
  subroutine DescribeDecomposition(DataMap,docfile)
  implicit none
  
  type( DataDecomposition )   :: DataMap
  character(len=80)           :: docfile
  integer                     :: myid
  integer                     :: numprocs
  integer                     :: ierr
  
  include 'mpif.h'
  
  call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) 
  call mpi_comm_size(MPI_COMM_WORLD,numprocs,ierr)  
  
  if(myid .ne. 0 ) return
  open(1,file=docfile,position='append')
  write(1,*) '============================================================'
  write(1,*) ' Data Decomposition ',trim(DataMap%label),' for : ',DataMap%numprocs,'   processors'
  write(1,*) '============================================================'

  write(1,*) ' Grid Level: ',DataMap%level
  write(1,*) ' Global number of data points in the '
  write(1,*) ' (x,y,z)<-->(1,2,3)<-->(i,j,k) directions:'
  write(1,*) '        nx = ',DataMap%nx
  write(1,*) '        ny = ',DataMap%ny
  write(1,*) '        nz = ',DataMap%nz
  write(1,*) ' Derivative types specified for each dimension: '
  write(1,*) '        in x ', trim(DataMap%DerivTypes(1))
  write(1,*) '        in y ', trim(DataMap%DerivTypes(2))
  write(1,*) '        in z ', trim(DataMap%DerivTypes(3))
  write(1,*) ' In {slab} format, each processor is allocated '
  write(1,*) ' nx*ny*locnz data values per variable, where'
  write(1,*) '        nx = ',DataMap%SlabVals(1,0)
  write(1,*) '        ny = ',DataMap%SlabVals(2,0)
  write(1,*) '        locnz = ',DataMap%SlabVals(3,0)
  if(DataMap%SlabVals(3,0) .ne. DataMap%SlabVals(3,DataMap%numprocs-1)) then
   write(1,*) ' except for processor ',DataMap%numprocs-1
   write(1,*) ' which is allocated 1 extra plane. '
   write(1,*) '       locnz = ',DataMap%SlabVals(3,DataMap%numprocs-1)
  endif
  write(1,*) '  '
  write(1,*) ' The local {slab} data arrays should be dimensioned                '
  write(1,*) '        nx = ',DataMap%SlabDims(1,0)
  write(1,*) '        ny = ',DataMap%SlabDims(2,0)
  write(1,*) '        locnz = ',DataMap%SlabDims(3,0)
  if(DataMap%SlabDims(3,0) .ne. DataMap%SlabDims(3,DataMap%numprocs-1)) then
   write(1,*) ' except for processor ',DataMap%numprocs-1
   write(1,*) ' which is allocated 1 extra plane. '
   write(1,*) '       locnz = ',DataMap%SlabVals(3,DataMap%numprocs-1)
   write(1,*) '  '
   write(1,*) ' myid = ',myid
   write(1,*) '   DataMap%SlabDims(1,myid) =     nx = ',DataMap%SlabDims(1,myid)
   write(1,*) '   DataMap%SlabDims(2,myid) =     ny = ',DataMap%SlabDims(2,myid)
   write(1,*) '   DataMap%SlabDims(3,myid) =     locnz = ',DataMap%SlabDims(3,myid)
   write(1,*) '   DataMap%locnz            =     locnz = ',DataMap%locnz
   write(1,*) '   DataMap%my_global_kvals  =             ',DataMap%my_global_kvals
   
   
  endif
  write(1,*) ' '
  write(1,*) '=========== '
  write(1,*) ' '
  write(1,*) ' In {column} format, measured in REAL words,'
  write(1,*) ' (even though the col arrays are sometimes complex),'
  write(1,*) ' each processor is allocated locnx*ny*nz '
  write(1,*) ' REAL data values per variable, where '
  write(1,*) '        locnx = ',DataMap%ColVals(1,0)
  write(1,*) '        ny = ',DataMap%ColVals(2,0)
  write(1,*) '        nz = ',DataMap%ColVals(3,0)
  if(DataMap%ColVals(1,0) .ne. DataMap%ColVals(1,DataMap%numprocs-1)) then
   write(1,*) ' except for processor ',DataMap%numprocs-1
   write(1,*) ' which is allocated 1 extra slot. '
   write(1,*) '       locnx = ',DataMap%ColVals(1,DataMap%numprocs-1)
  endif
  write(1,*) '  '
  write(1,*) ' The local {Col} data arrays should be dimensioned                '
  write(1,*) '        locnx = ',DataMap%ColVals(1,0)
  write(1,*) '        ny = ',DataMap%ColVals(2,0)
  write(1,*) '        nz = ',DataMap%ColVals(3,0)
  if(DataMap%ColVals(1,0) .ne. DataMap%ColVals(1,DataMap%numprocs-1)) then
   write(1,*) ' except for processor ',DataMap%numprocs-1
   write(1,*) ' which is allocated 1 extra plane. '
   write(1,*) '       locnx = ',DataMap%ColVals(1,DataMap%numprocs-1)
   write(1,*) '  '
   write(1,*) ' myid = ',myid
   write(1,*) '   DataMap%ColDims(1,myid) =     nx = ',DataMap%ColDims(1,myid)
   write(1,*) '   DataMap%ColDims(2,myid) =     ny = ',DataMap%ColDims(2,myid)
   write(1,*) '   DataMap%ColDims(3,myid) =     locnx = ',DataMap%ColDims(3,myid)
  endif
   write(1,*) '============================================================'
  write(1,*) '   '
  close(1) 
 end subroutine DescribeDecomposition
 

end Module class_DataDecomposition