/*=========================================================================

Program:   Visualization Toolkit
Module:    vtkPNetCDFPOPReader.cxx

Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#include "vtkPNetCDFPOPReader.h"
#include "vtkCallbackCommand.h"
#include "vtkDataArraySelection.h"
#include "vtkFloatArray.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkRectilinearGrid.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"

#include "vtkMPI.h"
#include "vtkMPIController.h"

#include "vtk_netcdf.h"
#include <string>
#include <vector>
#include <set>

#ifdef MPI_Comm
    #error MPI_Comm is #define'd somewhere!  That's BAD!  (Try checking netcdf.h.)
#endif

vtkStandardNewMacro(vtkPNetCDFPOPReader);

//============================================================================
#define CALL_NETCDF(call) \
{ \
  int errorcode = call; \
  if (errorcode != NC_NOERR) \
  { \
    vtkErrorMacro(<< "netCDF Error: " << nc_strerror(errorcode)); \
    return 0; \
  } \
}
//============================================================================


class vtkPNetCDFPOPReaderInternal
{
public:
  vtkSmartPointer<vtkDataArraySelection> VariableArraySelection;
  // a mapping from the list of all variables to the list of available
  // point-based variables
  std::vector<int> VariableMap;

  // Pointers to the buffers where we stored the values for each depth that we
  // read.  We can't delete the buffers until we know all the send requests that
  // used them (there may be more than one send for each buffer) have completed.
  std::vector <float *> SendBufs;

  // MPI ranks of the processes that will actually do the netCDF reads
  // NOTE: I'm currently getting these from an environment variable, but
  // I'm not certain that's the best way.  (If the var doesn't exist, I've
  // added what I hope are sane defaults.  See SetReaderRanks().)
  std::vector<int> ReaderRanks;

  // Memory to hold the extents for all processes (reader processes need this,
  // others can delete it after the Allgather operation, but it's still more
  // efficient to do an Allgather than to a bunch of individual Gathers.)
  int *AllExtents;

  // MPI_Request identifiers for all the MPI_Isend() calls.
  std::vector <MPI_Request> SendReqs;

  //////////////////////


  vtkPNetCDFPOPReaderInternal()
  {
      this->VariableArraySelection =
        vtkSmartPointer<vtkDataArraySelection>::New();
  }
};

//----------------------------------------------------------------------------
//set default values
vtkPNetCDFPOPReader::vtkPNetCDFPOPReader()
{
  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(1);
  this->FileName = NULL;
  this->OpenedFileName = NULL;
  this->Stride[0] = this->Stride[1] = this->Stride[2] = 1;
  this->SelectionObserver = vtkCallbackCommand::New();
  this->SelectionObserver->SetCallback
    (&vtkPNetCDFPOPReader::SelectionModifiedCallback);
  this->SelectionObserver->SetClientData(this);
  this->Internals = new vtkPNetCDFPOPReaderInternal;
  this->Internals->VariableArraySelection->AddObserver(
    vtkCommand::ModifiedEvent, this->SelectionObserver);
  this->Controller = NULL;
  this->SetController(vtkMPIController::SafeDownCast(
                        vtkMultiProcessController::GetGlobalController()));
  this->NCDFFD = -1;
}

//----------------------------------------------------------------------------
//delete filename and netcdf file descriptor
vtkPNetCDFPOPReader::~vtkPNetCDFPOPReader()
{
  this->SetController(NULL);
  this->SetFileName(0);
  if(this->OpenedFileName)
  {
    nc_close(this->NCDFFD);
  }
  this->SetOpenedFileName(0);
  if(this->SelectionObserver)
  {
    this->SelectionObserver->Delete();
    this->SelectionObserver = NULL;
  }
  delete this->Internals;
  this->Internals = NULL;
}

//----------------------------------------------------------------------------
void vtkPNetCDFPOPReader::PrintSelf(ostream &os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "FileName: "
     << (this->FileName ? this->FileName : "(NULL)") << endl;
  os << indent << "OpenedFileName: "
     << (this->OpenedFileName ? this->OpenedFileName : "(NULL)") << endl;
  os << indent << "Stride: {" << this->Stride[0] << ", "
     << this->Stride[1] << ", " << this->Stride[2] << ", "
     << "}" << endl;
  if(this->Controller)
  {
    os << indent << "Controller: " << this->Controller << endl;
  }
  else
  {
    os << indent << "Controller: (NULL)" << endl;
  }
  os << indent << "NCDFFD: " << this->NCDFFD << endl;

  this->Internals->VariableArraySelection->PrintSelf(os, indent.GetNextIndent());
}

//----------------------------------------------------------------------------
// RequestInformation supplies global meta information
// This should return the reality of what the reader is going to supply.
// This retrieve the extents for the rectilinear grid
// NC_MAX_VAR_DIMS comes from the nc library
int vtkPNetCDFPOPReader::RequestInformation(
    vtkInformation* vtkNotUsed(request),
    vtkInformationVector** vtkNotUsed(inputVector),
    vtkInformationVector* outputVector)
{
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  if(this->FileName == NULL)
  {
    vtkErrorMacro("FileName not set.");
    return 0;
  }

  // every rank that is a reader process needs to open the file
  if (this->IsReaderRank())
  {

    if(this->OpenedFileName)
    {
      nc_close(this->NCDFFD);
    }
    int retval = nc_open(this->FileName, NC_NOWRITE, &this->NCDFFD);//read file
    if (retval != NC_NOERR)//checks if read file error
    {
      vtkErrorMacro(<< "Can't read file " << nc_strerror(retval));
      this->SetOpenedFileName(0);
      return 0;
    }
    this->SetOpenedFileName(this->FileName);
  }
  // first reader reads the metadata and broadcasts it to everyone else
  int extent[6];
  char variableName[NC_MAX_NAME+1];
  if (this->IsFirstReaderRank())
  {
    // get number of variables from file
    int numberOfVariables;
    nc_inq_nvars(this->NCDFFD, &numberOfVariables);
    int dimidsp[NC_MAX_VAR_DIMS];
    int dataDimension;
    size_t dimensions[4]; //dimension value
    this->Internals->VariableMap.resize(numberOfVariables);
    int actualVariableCounter = 0;
    // For every variable in the file
    for(int i=0;i<numberOfVariables;i++)
    {
      this->Internals->VariableMap[i] = -1;
      //get number of dimensions
      CALL_NETCDF(nc_inq_varndims(this->NCDFFD, i, &dataDimension));
      //Variable Dimension ID's containing x,y,z coords for the rectilinear
      //grid spacing
      CALL_NETCDF(nc_inq_vardimid(this->NCDFFD, i, dimidsp));
      if(dataDimension == 3)
      {
        this->Internals->VariableMap[i] = actualVariableCounter++;
        //get variable name
        CALL_NETCDF(nc_inq_varname(this->NCDFFD, i, variableName));
        this->Internals->VariableArraySelection->AddArray(variableName);
        for(int m=0;m<dataDimension;m++)
        {
          CALL_NETCDF(nc_inq_dimlen(this->NCDFFD, dimidsp[m], dimensions+m));
          //acquire variable dimensions
        }
        extent[0] = extent[2] = extent[4] =0; //set extent
        extent[1] = static_cast<int>((dimensions[2] -1) / this->Stride[0]);
        extent[3] = static_cast<int>((dimensions[1] -1) / this->Stride[1]);
        extent[5] = static_cast<int>((dimensions[0] -1) / this->Stride[2]);
      }
    }

    // We've read in all the metadata.  Now broadcast it to the other ranks

    // There's probably only 1 variable name, but we'll allow for more just in case
    int numNames = this->Internals->VariableArraySelection->GetNumberOfArrays();
    this->Controller->Broadcast( &numNames, 1,
                                 this->Internals->ReaderRanks[0]);
    for (int i = 0; i < numNames; i++)
    {
      // I don't really like this extra strcpy, but I want to broadcast a string of
      // known, fixed size and the pointer retrurned by GetArrayName could point to
      // considerably less than NC_MAX_NAME chars, which means I'd risk a segfault if
      // I didn't do the copy.
      strcpy( variableName, this->Internals->VariableArraySelection->GetArrayName( i));
      this->Controller->Broadcast( variableName, NC_MAX_NAME+1,
                                   this->Internals->ReaderRanks[0]);
    }

    // Send out the variable map data
    int numVariables = static_cast<int>(this->Internals->VariableMap.size());
    this->Controller->Broadcast( &numVariables, 1,
                                 this->Internals->ReaderRanks[0]);
    this->Controller->Broadcast( &this->Internals->VariableMap[0], numVariables,
                                 this->Internals->ReaderRanks[0]);

    // send out the extents data
    this->Controller->Broadcast( extent, 6, this->Internals->ReaderRanks[0]);
  }
  else  // everyone else listens for the broadcasted metadata and fills in
        // Internals->VariableMap and extent[]
  {
    // Receive the variable name(s)  (probably only 1 name)
    int numNames = 0;
    this->Controller->Broadcast( &numNames, 1,
                                 this->Internals->ReaderRanks[0]);
    for (int i=0; i < numNames; i++)
    {
      this->Controller->Broadcast( variableName, NC_MAX_NAME+1,
                                   this->Internals->ReaderRanks[0]);
      this->Internals->VariableArraySelection->AddArray(variableName);
    }

    // Receive the variable map data
    int numVariables = 0;
    this->Controller->Broadcast( &numVariables, 1,
                                 this->Internals->ReaderRanks[0]);
    this->Internals->VariableMap.resize( numVariables, 0);
    this->Controller->Broadcast( &this->Internals->VariableMap[0], numVariables,
                                 this->Internals->ReaderRanks[0]);
    // Receive the extents data
    this->Controller->Broadcast( extent, 6, this->Internals->ReaderRanks[0]);
  }

  //fill in the extent information
  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
  outInfo->Set(CAN_PRODUCE_SUB_EXTENT(), 1);

  return 1;
}

//----------------------------------------------------------------------------
// Setting extents of the rectilinear grid
int vtkPNetCDFPOPReader::RequestData(vtkInformation* request,
    vtkInformationVector** vtkNotUsed(inputVector),
    vtkInformationVector* outputVector  )
{
  this->UpdateProgress(0);
  // the default implementation is to do what the old pipeline did find what
  // output is requesting the data, and pass that into ExecuteData
  // which output port did the request come from
  int outputPort = request->Get(vtkDemandDrivenPipeline::FROM_OUTPUT_PORT());
  // if output port is negative then that means this filter is calling the
  // update directly, in that case just assume port 0
  if (outputPort == -1)
  {
    outputPort = 0;
  }

  // get the data object
  vtkInformation *outInfo = outputVector->GetInformationObject(outputPort);
  vtkDataObject* output = outInfo->Get(vtkDataObject::DATA_OBJECT());

  int subext[6];
  //vtkInformation * outInfo = output->GetInformationObject(0);
  outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),subext);

  vtkRectilinearGrid *rgrid = vtkRectilinearGrid::SafeDownCast(output);
  rgrid->SetExtent(subext);

  size_t count[]= { static_cast<size_t>(subext[5]-subext[4]+1),
                    static_cast<size_t>(subext[3]-subext[2]+1),
                    static_cast<size_t>(subext[1]-subext[0]+1) };


  //initialize memory (raw data space, x y z axis space) and rectilinear grid
  bool firstPass = true;
  for(size_t i=0;i<this->Internals->VariableMap.size();i++)
  {
    if(this->Internals->VariableMap[i] != -1 &&
       this->Internals->VariableArraySelection->GetArraySetting(
         this->Internals->VariableMap[i]) != 0)
    {
      // varidp is probably i in which case nc_inq_varid isn't needed
      int varidp = VTK_INT_MIN;
      if (IsReaderRank())
      {
        nc_inq_varid(this->NCDFFD,
                     this->Internals->VariableArraySelection->GetArrayName(
                     this->Internals->VariableMap[i]), &varidp);
      }

      if(firstPass == true)
      {
        firstPass = false;
        // Get the latitude, longitude & depth values: the first reader process
        // does the read and broadcasts to everyone

        //set up start and stride values for this process's data (count is already set)
        size_t start[]= { static_cast<size_t>(subext[4]),
                          static_cast<size_t>(subext[2]),
                          static_cast<size_t>(subext[0]) };

        ptrdiff_t rStride[3] = { (ptrdiff_t)this->Stride[2],
                                 (ptrdiff_t)this->Stride[1],
                                 (ptrdiff_t)this->Stride[0] };

        int wholeExtent[6];
        outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),wholeExtent);
        // The Z & X dimensions have been swapped in wholeExtent (see the comments at the top
        // of the header file).  We want the arrays below to reflect the on-disk layout, which
        // is why the array indexes don't match up

        size_t wholeCount[3] = { static_cast<size_t>(wholeExtent[5] - wholeExtent[4] + 1),
                                 static_cast<size_t>(wholeExtent[3] - wholeExtent[2] + 1),
                                 static_cast<size_t>(wholeExtent[1] - wholeExtent[0] + 1) };

        size_t wholeStart[3]= { static_cast<size_t>(wholeExtent[4]*this->Stride[2]),
                                static_cast<size_t>(wholeExtent[2]*this->Stride[1]),
                                static_cast<size_t>(wholeExtent[0]*this->Stride[0]) };

        std::vector<float> buffer(wholeCount[0] + wholeCount[1] + wholeCount[2]);
        if (this->IsFirstReaderRank())
        {
          int dimidsp[3];
          nc_inq_vardimid(this->NCDFFD, varidp, dimidsp);
          nc_get_vars_float(this->NCDFFD, dimidsp[0], wholeStart, wholeCount,
                            rStride, &buffer[0]);
          nc_get_vars_float(this->NCDFFD, dimidsp[1], wholeStart+1, wholeCount+1,
                            rStride+1, &buffer[wholeCount[0]]);
          nc_get_vars_float(this->NCDFFD, dimidsp[2], wholeStart+2, wholeCount+2,
                            rStride+2, &buffer[wholeCount[0] + wholeCount[1]]);
        }

        this->Controller->Broadcast( &buffer[0], (wholeCount[0] + wholeCount[1] + wholeCount[2]),
                                     this->Internals->ReaderRanks[0]);

        // Extract the values we need out of buffer and store them in the x, y & z buffers
        float* x = new float[count[0]];
        float* y = new float[count[1]];
        float* z = new float[count[2]];
        memcpy( x, &buffer[0 + start[0]],                            count[0] * sizeof( float));
        memcpy( y, &buffer[wholeCount[0] + start[1]],                count[1] * sizeof( float));
        memcpy( z, &buffer[wholeCount[0] + wholeCount[1]+ start[2]], count[2] * sizeof( float));

        // Note the axis swap:  xcoords gets the z data and zcoords gets the x data
        vtkFloatArray *xCoords = vtkFloatArray::New();
        xCoords->SetArray(z, count[2], 0, 1);
        vtkFloatArray *yCoords = vtkFloatArray::New();
        yCoords->SetArray(y, count[1], 0, 1);
        for (unsigned int q=0; q<count[0]; q++)
        {
          x[q] = -x[q];
        }
        vtkFloatArray *zCoords = vtkFloatArray::New();
        zCoords->SetArray(x, count[0], 0, 1);
        rgrid->SetXCoordinates(xCoords);
        rgrid->SetYCoordinates(yCoords);
        rgrid->SetZCoordinates(zCoords);
        xCoords->Delete();
        yCoords->Delete();
        zCoords->Delete();
      }
      //create vtkFloatArray and get the scalars into it
      vtkFloatArray *scalars = vtkFloatArray::New();
      vtkIdType numberOfTuples = (count[0])*(count[1])*(count[2]);
      float* data = new float[numberOfTuples];


      // Do a gather of all process' sub-extents so that the reader processes will know
      // who needs what data.  (An AllGather() operation is somewhat wasteful, because
      // even processes that aren't readers will still receive the data, but it's still
      // more efficient than having everyone send extents to each individual reader process.)
      int mpiNumProcs = this->Controller->GetNumberOfProcesses();
      this->Internals->AllExtents = new int[ 6 * mpiNumProcs];
      this->Controller->AllGather(subext, this->Internals->AllExtents, 6);

      // First, post all the receive requests
      std::vector <MPI_Request> recvReqs;  // should be 1 recv for each depth

      // Number of values stored for each depth
      //unsigned long oneDepthSize = (subext[3]-subext[2]+1) * (subext[5]-subext[4]+1);
      unsigned long oneDepthSize = static_cast<unsigned long>(count[1]*count[2]);  // should be the same value as the line above...

      for (int curDepth = subext[4]; curDepth <= subext[5]; curDepth++)
      {
        float *depthStart = data + ((curDepth-subext[4])*oneDepthSize);  // Where this data should start
        int sourceRank =  ReaderForDepth( curDepth);
        //      cerr << "Rank " << mpiRank << ": Expecting depth " << curDepth << " from rank " << sourceRank << endl;
        // We could probably use vtkMPIController::NoBlockReceive(), but but I'm not sure
        // how it will behave since we're calling MPI_Isend() directly.  So, I'm going
        // to call MPI_Irecv directly, too...

        MPI_Comm *comm = ((vtkMPICommunicator *)this->Controller->GetCommunicator())->GetMPIComm()->GetHandle();
        MPI_Request recvReq;
        MPI_Irecv( depthStart, oneDepthSize, MPI_FLOAT, sourceRank, curDepth, *comm, &recvReq);
        recvReqs.push_back( recvReq);
      }


      if (IsReaderRank())
      {
        // Reads part of the netCDF file and sends subarrays out to all the ranks
        ReadAndSend( outInfo, varidp);
      }

      delete[] this->Internals->AllExtents;


      // Wait for all the sends to complete
      if (this->Internals->SendReqs.size() > 0)
      {
        MPI_Waitall( static_cast<int>(this->Internals->SendReqs.size()),
                     &this->Internals->SendReqs[0], MPI_STATUSES_IGNORE);

        // Now that all the sends are complete, it's safe to free the buffers
        for (size_t j= 0; j < this->Internals->SendBufs.size(); j++)
        {
          delete[] this->Internals->SendBufs[j];
        }
        this->Internals->SendBufs.clear();
        this->Internals->SendReqs.clear();
      }

      MPI_Waitall( static_cast<int>(recvReqs.size()), &recvReqs[0], MPI_STATUSES_IGNORE);
      recvReqs.clear();


      scalars->SetArray(data, numberOfTuples, 0, 1);
      //set list of variables to display data on rectilinear grid
      scalars->SetName(this->Internals->VariableArraySelection->GetArrayName(
                         this->Internals->VariableMap[i]));
      rgrid->GetPointData()->AddArray(scalars);
      scalars->Delete();
    }
    this->UpdateProgress((i+1.0)/this->Internals->VariableMap.size());
  }
  return 1;
}

//----------------------------------------------------------------------------
//following 5 functions are used for paraview user interface
void vtkPNetCDFPOPReader::SelectionModifiedCallback(vtkObject*, unsigned long,
                                                   void* clientdata, void*)
{
  static_cast<vtkPNetCDFPOPReader*>(clientdata)->Modified();
}

//-----------------------------------------------------------------------------
int vtkPNetCDFPOPReader::GetNumberOfVariableArrays()
{
  return this->Internals->VariableArraySelection->GetNumberOfArrays();
}

//-----------------------------------------------------------------------------
const char* vtkPNetCDFPOPReader::GetVariableArrayName(int index)
{
  if(index < 0 || index >= this->GetNumberOfVariableArrays())
  {
    return NULL;
  }
  return this->Internals->VariableArraySelection->GetArrayName(index);
}

//-----------------------------------------------------------------------------
int vtkPNetCDFPOPReader::GetVariableArrayStatus(const char* name)
{
  return this->Internals->VariableArraySelection->ArrayIsEnabled(name);
}

//-----------------------------------------------------------------------------
void vtkPNetCDFPOPReader::SetVariableArrayStatus(const char* name, int status)
{
  vtkDebugMacro("Set cell array \"" << name << "\" status to: " << status);
  if(this->Internals->VariableArraySelection->ArrayExists(name) == 0)
  {
    vtkErrorMacro(<< name << " is not available in the file.");
    return;
  }
  int enabled = this->Internals->VariableArraySelection->ArrayIsEnabled(name);
  if(status != 0 && enabled == 0)
  {
    this->Internals->VariableArraySelection->EnableArray(name);
    this->Modified();
  }
  else if(status == 0 && enabled != 0)
  {
    this->Internals->VariableArraySelection->DisableArray(name);
    this->Modified();
  }
}

namespace
{
  // shortcut used down in ReadAndSend()
  void swap( int &A, int &B) { int temp = B; B=A; A=temp;}
}

//----------------------------------------------------------------------------
// Helper function for RequestData.  Reads one or more depth arrays from the
// netCDF file and sends sub-arrays out to all ranks that need that data
int vtkPNetCDFPOPReader::ReadAndSend( vtkInformation *outInfo, int varID)
{
  int wholeExtent[6];
  outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),wholeExtent);
  // Z & X dimensions have been swapped (see the comments at the top of the header
  // file).  All arrays in this function, however, should reflect the on-disk
  // layout, so we'll 'un-swap' the X & Z extents.
  swap( wholeExtent[0], wholeExtent[4]);
  swap( wholeExtent[1], wholeExtent[5]);

  // this->Stride is also in the in-memory layout
  ptrdiff_t rStride[3] = { this->Stride[2], this->Stride[1], this->Stride[0] };

  int rank = this->Controller->GetLocalProcessId();

  // We read one depth at a time, skipping over the depths that other reader processes will read
  for (int curDepth = wholeExtent[0]; curDepth <= wholeExtent[1]; curDepth++)
  {
    if (ReaderForDepth( curDepth) == rank)
    {
      size_t start[3] = { static_cast<size_t>(curDepth*rStride[0]),
                          static_cast<size_t>(wholeExtent[2]),
                          static_cast<size_t>(wholeExtent[4]) };
      size_t count[3] = { 1,
                          static_cast<size_t>(wholeExtent[3] - wholeExtent[2] + 1),
                          static_cast<size_t>(wholeExtent[5] - wholeExtent[4] + 1) };

      float* buffer = new float[count[1]*count[2]];
      this->Internals->SendBufs.push_back( buffer );

      int ncErr = nc_get_vars_float( this->NCDFFD, varID, start, count, rStride, buffer);
      if (ncErr != NC_NOERR)
      {
        cerr << "!!!nc_get_vars_float() returned error code " << ncErr << endl;
      }

      // Create sub arrays and send to all processes
      for (int destRank = 0; destRank < this->Controller->GetNumberOfProcesses(); destRank++)
      {
        int destExtent[6];
        memcpy( destExtent, &this->Internals->AllExtents[destRank*6], 6*sizeof( int));
        // Note that AllExtents is also in in-memory layout order, so we need to swap
        // the X & Z values
        swap( destExtent[0], destExtent[4]);
        swap( destExtent[1], destExtent[5]);

        // Verify that destRank does, in fact, receive an extent at this depth
        // TODO: Don't use MPI to send data to ourself.....
        if (curDepth >= destExtent[0] && curDepth <= destExtent[1])
        {
          int subarray_sizes[2] =    { (wholeExtent[3] - wholeExtent[2]+1), (wholeExtent[5] - wholeExtent[4]+1), };
          int subarray_subsizes[2] = { (destExtent[3] -  destExtent[2]+1),  (destExtent[5]  - destExtent[4]+1) };
          int subarray_starts[2] = { destExtent[2], destExtent[4] };
          MPI_Datatype subArrayType;
          MPI_Request sendReq;
          MPI_Type_create_subarray( 2, subarray_sizes, subarray_subsizes, subarray_starts, MPI_ORDER_C, MPI_FLOAT, &subArrayType);
          MPI_Type_commit( &subArrayType);
          // cerr << "Rank " << mpiRank << ": Sending depth " << curDepth << " to rank " << destRank << endl;

          // vtkMPICommunicator can't handle arbitrary types, so we'll have to do it ourselves
          MPI_Comm *comm = ((vtkMPICommunicator *)this->Controller->GetCommunicator())->GetMPIComm()->GetHandle();
          MPI_Isend( buffer, 1, subArrayType, destRank, curDepth, *comm, &sendReq);  // using the depth value as the tag
          this->Internals->SendReqs.push_back( sendReq);
          MPI_Type_free( &subArrayType);
        }
      }

      // Check to see if any of the previous sends have completed
      // (This helps to keep the number of 'in-flight' sends to a minimum.)
      if (this->Internals->SendReqs.size())
      {
        int foundOne = 1;
        int reqIndex = 0;
        MPI_Status status; // we never actually use this, but MPI_Testany() requires it

        // MPI_Testany() will deallocate the request, so we don't need to do anything
        // except call it in a loop.
        do
        {
          MPI_Testany( static_cast<int>(this->Internals->SendReqs.size()),
                       &this->Internals->SendReqs[0],
                       &reqIndex, &foundOne, &status);
        } while (foundOne && reqIndex != MPI_UNDEFINED);
      }
    }
  }

  return 1;
}


//----------------------------------------------------------------------------
// Sets the ranks that will actually open and read the netcdf file.  If the
// ranks pointer is null, picks a default set.  The end result is an updated
// readerRanks vector in the Internals object.
void vtkPNetCDFPOPReader::SetReaderRanks(vtkIdList* ranks)
{
  std::set<int>readerRanks;  // temporary holder for the values.
  // We use a set so that ranks are automatically ordered and duplicates
  // are skipped.  Once the set is created, it's copied over to a vector
  // because I want to be able to use operator[] to access the values.

  int numProcs = this->Controller->GetNumberOfProcesses();

  if(ranks)
  {
    for(vtkIdType i=0;i<ranks->GetNumberOfIds();i++)
    {
      vtkIdType rank = ranks->GetId(i);
      if(rank >= 0 && rank < numProcs)
      {
        readerRanks.insert(rank);
      }
    }
  }

  if (readerRanks.empty())  // Either nobody set the env var or it had bogus values
                            //  in it.  Try to pick a reasonable default.
  {
    // This is somewhat arbritrary:  below 24 processes, we'll use 4 readers.
    // More than 24 processes, we'll use 8.  (>24 processes, I'm assuming we're
    // running on Jaguar where 2 readers per OSS seems to work well).
    // All readers will be evenly spread across the range of processes that
    // are working on this file

    int numReaders = numProcs < 24 ? 4 : 8;

    if (numProcs < numReaders)  // sanity check
    {
      numReaders = numProcs;
    }
    int proc = 0;
    while (proc < numProcs)
    {
      readerRanks.insert( proc);
      proc += (numProcs / numReaders);
    }
  }

  // Copy the contents of the set into its permanent location...
  this->Internals->ReaderRanks.clear();
  std::set<int>::iterator it = readerRanks.begin();
  while (it != readerRanks.end())
  {
    this->Internals->ReaderRanks.push_back( *it);
    it++;
  }
}

//----------------------------------------------------------------------------
// Returns the rank (relative to our controller) for the process that will
// read the specified depth
int vtkPNetCDFPOPReader::ReaderForDepth( unsigned depth)
{
  // NOTE: This is a very simple algorithm - each rank in readerRanks will
  // read single depth in a round-robbin fashion.  There might be a more
  // efficient way to do this...
  size_t numReaders = this->Internals->ReaderRanks.size();
  return this->Internals->ReaderRanks[(depth % numReaders)];
}

//----------------------------------------------------------------------------
// Returns true if the calling process should read data from the netCDF file
bool vtkPNetCDFPOPReader::IsReaderRank()
{
  int rank = this->Controller->GetLocalProcessId();
  for (size_t i = 0; i < this->Internals->ReaderRanks.size(); i++)
  {
    int ii = static_cast<int>(i);
    if (rank == this->Internals->ReaderRanks[ii])
    {
      return true;
    }
  }
  return false;
}

//----------------------------------------------------------------------------
// Similar to above, but returns true only if the calling process is the first
// rank in the readerRanks vector.  (This function exists because much of the
// file metadata is read by a single rank and broadcast to all the others.)
bool vtkPNetCDFPOPReader::IsFirstReaderRank()
{
  if (this->Internals->ReaderRanks.size() == 0)
  {
    return false; // sanity check
  }

  int rank = this->Controller->GetLocalProcessId();
  return (rank == this->Internals->ReaderRanks[0]);
}
//----------------------------------------------------------------------------
//
void vtkPNetCDFPOPReader::SetController(vtkMPIController *controller)
{
  if(this->Controller != controller)
  {
    this->Controller = controller;
    if (this->Controller != NULL)
    {
      this->SetReaderRanks(NULL);
    }
  }
}