PLAG_RFLO_ModStats.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Collection of routines for particle statistics on Eulerian grid.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: PLAG_RFLO_ModStats.F90,v 1.9 2009/03/02 00:19:36 mtcampbe Exp $
!
! Copyright: (c) 2004 by the University of Illinois
!
! ******************************************************************************

MODULE plag_rflo_modstats

  USE modparameters
  USE moddatatypes  
  USE modglobal, ONLY: t_global
  USE modpartlag, ONLY: t_buffer_plag
  USE modbndpatch, ONLY: t_patch
  USE modgrid, ONLY: t_grid
  USE modpartlag, ONLY: t_plag 
  USE moddatastruct, ONLY: t_region
  USE moderror
  USE modmpi
  USE plag_modparameters

#include "Indexing.h" 
  USE modindexing,   ONLY: indijkmap
  USE modinterfaces, ONLY: rflo_getcelloffset,     &
                           rflo_getdimensdummy,    &
                           rflo_getpatchindices,   &
                           rflo_getpatchdirection, &
                           rflo_getpatchmapping,   &
                           rflo_readdatafileint,   &
                           rflo_readdatafilereal,  &
                           rflo_writedatafileint,  &
                           rflo_writedatafilereal

  IMPLICIT NONE
  
! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

! ==============================================================================
! Private data
! ==============================================================================

  CHARACTER(CHRLEN), PRIVATE :: & 
    RCSIdentString = '$RCSfile: PLAG_RFLO_ModStats.F90,v $ $Revision: 1.9 $'        

! ==============================================================================
! Public functions
! ==============================================================================
  
  PUBLIC :: plag_rflo_createstatbuff,      &
            plag_rflo_commstatbuffwrapper, &
            plag_rflo_readstat,            &
            plag_rflo_writestat

! ==============================================================================
! Private functions
! ==============================================================================

! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  
  
  
  
  
  
  
  
  
  

  





!******************************************************************************
!
! Purpose: Clear communication requests.
!
! Description: none.
!
! Input: regions = data of all regions
!        iReg    = index of current region.
!
! Output: none.
!
! Notes: none.
!
!******************************************************************************

  SUBROUTINE plag_rflo_clearreqstatbuff( regions,iReg )

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: regions(:)
    INTEGER :: ireg

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: bctype,ilev,ipatch,iregsrc,irequeststat,npatches 
#ifdef MPI
    INTEGER :: status(mpi_status_size)
#endif
    LOGICAL :: dowait

    TYPE(t_patch),    POINTER :: ppatch
    TYPE(t_plag),     POINTER :: pplag
    TYPE(t_region),   POINTER :: pregion
    TYPE(t_global),   POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(ireg)%global

    CALL registerfunction( global,'PLAG_RFLO_ClearReqStatBuff',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Set dimensions, variables, pointers
! ******************************************************************************    

#ifdef MPI
    ilev     = regions(ireg)%currLevel
    npatches = regions(ireg)%nPatches

    pplag => regions(ireg)%levels(ilev)%plag 

! ****************************************************************************** 
!   Wait for patch data being received by other processors
! ******************************************************************************

    DO ipatch=1,npatches
      ppatch => regions(ireg)%levels(ilev)%patches(ipatch)

      bctype  = ppatch%bcType
      iregsrc = ppatch%srcRegion
      irequeststat = ppatch%bufferPlag%iRequestStat

      dowait = ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
                (bctype>=bc_regionint  .AND. bctype<=bc_regionint +bc_range) .OR. &
                (bctype>=bc_regnonconf .AND. bctype<=bc_regnonconf+bc_range) .OR. &
                (bctype>=bc_tra_peri   .AND. bctype<=bc_tra_peri  +bc_range) .OR. &
                (bctype>=bc_rot_peri   .AND. bctype<=bc_rot_peri  +bc_range))

      IF ( iregsrc > 0 ) THEN
        IF ( (dowait .EQV. .true.) .AND.                  &
             (regions(iregsrc)%procid /= global%myProcid) ) THEN
          CALL mpi_wait( pplag%requestsStat(irequeststat),status,global%mpierr )
          IF ( global%mpierr /= err_none ) &
            CALL errorstop( global,err_mpi_trouble,__line__ )
        ENDIF ! doWait 
      ENDIF ! iRegSrc

    ENDDO   ! iPatch
#endif

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_clearreqstatbuff






! ******************************************************************************
!
! Purpose: Create buffer arrays for statistics.
!
! Description: None.
!
! Input: regions = dimensions of all regions
!        iReg    = current region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_rflo_createstatbuff(regions,iReg)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: regions(:)
    INTEGER :: ireg 

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: bctype,errorflag,ilev,ipatch,iregsrc,n1,n2,n1src,n2src,ndc,&
               ndcsrc,ndim,ndimsrc,neqs,neqssrc,nev,ntav

    TYPE(t_patch),       POINTER :: ppatch
    TYPE(t_buffer_plag), POINTER :: pbuffplag
    TYPE(t_plag),        POINTER :: pplag
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(ireg)%global

    CALL registerfunction( global,'PLAG_RFLO_CreateStatBuff',&
  'PLAG_RFLO_ModStats.F90' )

    IF ( global%myProcid == masterproc .AND. & 
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,3X,A)') solver_name,&
                               'Allocating Statistics Data Buffers for PLAG...'
    END IF ! global%verbLevel
    
    ntav = global%plagNStat

! ******************************************************************************
!   Allocate buffer data arrays for statistics
! ******************************************************************************  

    DO ilev=1,regions(ireg)%nGridLevels
      pplag => regions(ireg)%levels(ilev)%plag
      pplag%nRequestsStat = 0
   
      nev = pplag%nEv
      
      DO ipatch=1,regions(ireg)%nPatches
  
        ppatch => regions(ireg)%levels(ilev)%patches(ipatch)
        pbuffplag => ppatch%bufferPlag
      
        bctype = ppatch%bcType

        IF ( (bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
             (bctype>=bc_tra_peri   .AND. bctype<=bc_tra_peri  +bc_range) .OR. &
             (bctype>=bc_rot_peri   .AND. bctype<=bc_rot_peri  +bc_range)) THEN
          iregsrc = ppatch%srcRegion
        
          IF ( regions(iregsrc)%procid /= global%myProcid ) THEN  ! other processor
            n1      = abs(ppatch%l1end   -ppatch%l1beg   ) + 2    ! large enough
            n2      = abs(ppatch%l2end   -ppatch%l2beg   ) + 2    ! for NODES!
            n1src   = abs(ppatch%srcL1end-ppatch%srcL1beg) + 2
            n2src   = abs(ppatch%srcL2end-ppatch%srcL2beg) + 2
            neqs    = ntav +nev     
            neqssrc = ntav +nev
            ndc     = regions(ireg   )%nDumCells
            ndcsrc  = regions(iregsrc)%nDumCells
            ndim    = n1*n2*neqs*ndc
            ndimsrc = n1src*n2src*neqssrc*ndcsrc

            ALLOCATE( pbuffplag%sendBuffStat(ndimsrc),stat=errorflag )
            global%error = errorflag
            IF (global%error /= err_none) THEN
              CALL errorstop( global, err_allocate,__line__,'pBuffPlag%sendBuffStat' ) 
            END IF ! global%error

            ALLOCATE( pbuffplag%recvBuffStat(ndim   ),stat=errorflag )
            global%error = errorflag
            IF (global%error /= err_none) THEN
              CALL errorstop( global, err_allocate,__line__,'pBuffPlag%recvBuffStat' ) 
            END IF ! global%error

            pbuffplag%nSendBuffStat = ndimsrc
            pbuffplag%nRecvBuffStat = ndim
            pplag%nRequestsStat     = pplag%nRequestsStat + 1
            pbuffplag%iRequestStat  = pplag%nRequestsStat

          ENDIF ! regions(iRegSrc)%procid 

        ELSE IF ( (bctype>=bc_regionint  .AND. bctype<=bc_regionint +bc_range) .OR. &
                  (bctype>=bc_regnonconf .AND. bctype<=bc_regnonconf+bc_range)) THEN
          CALL errorstop( global,err_unknown_bc,__line__ )  ! #### TEMPORARY ####

        ELSE 
          nullify(pbuffplag%sendBuffStat)
          nullify(pbuffplag%recvBuffStat)
        ENDIF  ! bcType

      ENDDO    ! iPatch

! ==============================================================================
!     Allocate array for send requests
! ==============================================================================

      ALLOCATE( pplag%requestsStat(pplag%nRequestsStat),stat=errorflag )
      global%error = errorflag
      IF (global%error /= err_none) THEN
        CALL errorstop( global, err_allocate,__line__,'pPlag%requests' ) 
      END IF ! global%error

    ENDDO      ! iLev

! ******************************************************************************
!   Finalize
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. & 
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,3X,A)') solver_name,&
                           'Allocating Statistics Data Buffers for PLAG done...'
    END IF ! global%verbLevel

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_createstatbuff






  
! ******************************************************************************
!
! Purpose: Communicate buffer arrays for statistics.
!
! Description: None.
!
! Input: regions = dimensions of all regions
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_rflo_commstatbuffwrapper(regions)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: regions(:)

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: bctype,ilev,ipatch,ipatchsrc,ireg,iregsrc,npatches
  
    TYPE(t_patch),       POINTER :: ppatch,ppatchsrc
    TYPE(t_region),      POINTER :: pregion,pregionsrc
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(1)%global

    CALL registerfunction( global,'PLAG_RFLO_CommStatBuffWrapper',&
  'PLAG_RFLO_ModStats.F90' )

    IF (.NOT. global%plagUsed) goto 999

! ******************************************************************************
!   Copy statistics for regions on the same processor
! ******************************************************************************  

    DO ireg=1,global%nRegions
      pregion => regions(ireg)
      
      IF ( pregion%procid==global%myProcid .AND. &   ! region active and
           pregion%active==active ) THEN             ! on my processor

! ==============================================================================
!        Set dimensions and pointers
! ==============================================================================

        ilev     = pregion%currLevel
        npatches = pregion%nPatches


! ==============================================================================         
!        Loop over patches 
! ==============================================================================
 
        DO ipatch=1,npatches
          ppatch => pregion%levels(ilev)%patches(ipatch)

          bctype    = ppatch%bcType
          iregsrc   = ppatch%srcRegion
          ipatchsrc = ppatch%srcPatch

          SELECT CASE (bctype)

! ------------------------------------------------------------------------------
!           Conforming region interface
! ------------------------------------------------------------------------------

            CASE( bc_regionconf:bc_regionconf+bc_range )
              pregionsrc => regions(iregsrc)
              ppatchsrc  => pregionsrc%levels(ilev)%patches(ipatchsrc)

              IF ( regions(iregsrc)%procid == global%myProcid ) THEN
                CALL plag_rflo_copystatbuff( pregion,pregionsrc, &
                                             ppatch,ppatchsrc    )
              ENDIF ! regions(iRegSrc)%procid

! ------------------------------------------------------------------------------
!           Non-conforming region interface (integer)
! ------------------------------------------------------------------------------

            CASE( bc_regionint:bc_regionint+bc_range )
              CALL errorstop( global,err_unknown_bc,__line__ )

! ------------------------------------------------------------------------------
!           Non-conforming region interface (irregular)
! ------------------------------------------------------------------------------
             
            CASE( bc_regnonconf:bc_regnonconf+bc_range)
              CALL errorstop( global,err_unknown_bc,__line__ )         
          
          END select! bcType 

        END DO ! iPatch
      END IF ! pRegion%procid
    END DO ! iReg      

! ******************************************************************************
!   Communicate buffer for off-processor regions
! ******************************************************************************  

! ==============================================================================
!  Send buffer data
! ==============================================================================

#ifdef MPI
    DO ireg = 1, global%nRegions
      IF ( regions(ireg)%procid==global%myProcid .AND. &   ! region active and
           regions(ireg)%active==active ) THEN             ! on my processor        
        CALL plag_rflo_sendstatbuffwrapper( regions, ireg )    
      ENDIF       ! regions    
    ENDDO         ! iReg

! ==============================================================================                     
!  Receive buffer data 
! ==============================================================================
 
    DO ireg = 1, global%nRegions
      IF ( regions(ireg)%procid==global%myProcid .AND. &   ! region active and
           regions(ireg)%active==active ) THEN             ! on my processor        
        CALL plag_rflo_recvstatbuffwrapper( regions, ireg )
      ENDIF ! regions
    ENDDO ! iReg 

! ==============================================================================
! wait for data being received by other processors ----------------------------
! ==============================================================================

    DO ireg = 1, global%nRegions
      IF ( regions(ireg)%procid==global%myProcid .AND. &   ! region active and
           regions(ireg)%active==active ) THEN             ! on my processor        
        CALL plag_rflo_clearreqstatbuff( regions, ireg )
      ENDIF ! regions  
    ENDDO ! iReg
#endif

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

999  CONTINUE
    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_commstatbuffwrapper






  
! ******************************************************************************
!
! Purpose: Copy buffer statistics arrays for on-processor regions.
!
! Description: None.
!
! Input: pRegion    = current region
!        pRegionSrc = source region
!        pPatch     = current patch
!        pPatchSrc  = source patch
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_rflo_copystatbuff( pRegion,pRegionSrc,pPatch,pPatchSrc )

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_patch),   POINTER :: ppatch,ppatchsrc
    TYPE(t_region),  POINTER :: pregion,pregionsrc

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: idum,iv,i,j,k,ii,jj,kk
    INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, idir, jdir, kdir, &
               icoff, ijcoff, ijkd, ilev
    INTEGER :: ibegsrc, iendsrc, jbegsrc, jendsrc, kbegsrc, kendsrc, &
               idirsrc, jdirsrc, kdirsrc, icoffsrc, ijcoffsrc, ijkcsrc
    INTEGER :: lb, lbs, l1srcdir, l2srcdir, mapmat(3,4)
    INTEGER :: ivevbeg,ivevend,ivtavbeg,ivtavend

    LOGICAL :: align
    
    REAL(KIND=RFREAL), DIMENSION(:,:), POINTER :: ev,evsrc,tav,tavsrc
  
    TYPE(t_plag),        POINTER :: pplag
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global

    CALL registerfunction( global,'PLAG_RFLO_CopyStatBuff',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Check if the source region is active
! ******************************************************************************
    
    IF ( pregionsrc%active == off ) THEN
      CALL errorstop( global,err_srcregion_off,__line__ )
    ENDIF

! ******************************************************************************
!   Set dimensions and pointers
! ******************************************************************************

    ilev = pregion%currLevel

    CALL rflo_getpatchindices( pregion,ppatch,ilev,ibeg,iend, &
                               jbeg,jend,kbeg,kend )
    CALL rflo_getpatchindices( pregionsrc,ppatchsrc,ilev,ibegsrc,iendsrc, &
                               jbegsrc,jendsrc,kbegsrc,kendsrc )
    CALL rflo_getpatchdirection( ppatch   ,idir   ,jdir,   kdir    )
    CALL rflo_getpatchdirection( ppatchsrc,idirsrc,jdirsrc,kdirsrc )
    CALL rflo_getcelloffset( pregion   ,ilev,icoff   ,ijcoff    )
    CALL rflo_getcelloffset( pregionsrc,ilev,icoffsrc,ijcoffsrc )
    
    ev    => pregion%levels(ilev)%plag%ev
    evsrc => pregionsrc%levels(ilev)%plag%ev
    
    ivevbeg = 1
    ivevend = SIZE(ev,dim=1)
        
    tav    => pregion%levels(ilev)%plag%tav
    tavsrc => pregionsrc%levels(ilev)%plag%tav
    
    ivtavbeg = 1
    ivtavend = SIZE(tav,dim=1)
    
! ******************************************************************************
!   Mapping between patches
! ******************************************************************************

                                           l1srcdir =  1
    IF (ppatch%srcL1beg > ppatch%srcL1end) l1srcdir = -1
                                           l2srcdir =  1
    IF (ppatch%srcL2beg > ppatch%srcL2end) l2srcdir = -1

    lb    = ppatch%lbound
    lbs   = ppatch%srcLbound
    align = ppatch%align

    CALL rflo_getpatchmapping( lb,lbs,l1srcdir,l2srcdir,align, &
                               idir,jdir,kdir,idirsrc,jdirsrc,kdirsrc, &
                               ibeg,iend,jbeg,jend,kbeg,kend, &
                               ibegsrc,iendsrc,jbegsrc,jendsrc,kbegsrc,kendsrc, &
                               mapmat )

! ******************************************************************************
!   Loop over dummy nodes of current patch
! ******************************************************************************
    
    DO idum=1,pregion%nDumCells
    DO k=kbeg,kend
    DO j=jbeg,jend
    DO i=ibeg,iend
      ii      = i - idum*idir
      jj      = j - idum*jdir
      kk      = k - idum*kdir
      ijkd    = indijk(ii,jj,kk,icoff,ijcoff)
      ijkcsrc = indijkmap(ii,jj,kk,mapmat,icoffsrc,ijcoffsrc)

      DO iv=ivevbeg,ivevend
         ev(iv,ijkd) =  evsrc(iv,ijkcsrc)
      END DO ! iVar

      DO iv=ivtavbeg,ivtavend
        tav(iv,ijkd) = tavsrc(iv,ijkcsrc)
      END DO ! iVar      
    ENDDO  ! i
    ENDDO  ! j
    ENDDO  ! k
    ENDDO  ! idum

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_copystatbuff






  

! ******************************************************************************
!
! Purpose: Read in time averaged statistics of the Lagrangian particles 
!          on Eulerian grid.
!
! Description: the following solution formats are supported:
!              - RocfloMP ASCII
!              - RocfloMP binary
!
! Input: regions = dimensions of all regions
!
! Output: region%levels%plag%tav = time avg Lagrangian particle variables
!         global%integrTime      = integrated averaging time
!
! Notes: time averaged solution is read in only for the current grid level;
!        it is also read in for all dummy cells
!
! ******************************************************************************

  SUBROUTINE plag_rflo_readstat(regions)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: regions(:)  

! ==============================================================================
!   Locals
! ==============================================================================

  CHARACTER(2*CHRLEN+17) :: fname
  CHARACTER(CHRLEN)      :: msg

  INTEGER :: ireg, i, j, k, l, n, ind

#ifdef MPI
  INTEGER :: status(mpi_status_size)
#endif
  INTEGER :: ilev,iregfile,ipc,jpc,kpc,ndumcells,ndim,ioff,ijoff,ijk
  INTEGER :: idcbeg,jdcbeg,kdcbeg,idcend,jdcend,kdcend,ijkbeg,ijkend
  INTEGER :: errorflag,ntav,ntavvar

  INTEGER, ALLOCATABLE, DIMENSION (:,:) :: ivar,jvar,plagvarid

  REAL(RFREAL), POINTER, DIMENSION(:,:)     :: tav
  REAL(RFREAL), ALLOCATABLE, DIMENSION(:,:) :: rvar, tavfile

  TYPE(t_global), POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(1)%global

    CALL registerfunction( global,'PLAG_RFLO_ReadStat',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Allocate temporary data arrays 
! ******************************************************************************

    ALLOCATE( ivar(5,1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'ivar' )
    
    ALLOCATE( rvar(2,1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'rvar' )
    
    ALLOCATE( jvar(global%plagNStat+1,1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'jvar' ) 
    
    ALLOCATE( plagvarid(2,global%plagNStat+1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'plagVarId' ) 

! ******************************************************************************
!   Open statistics file (only master proc.)
! ******************************************************************************

    IF (global%myProcid == masterproc) THEN

      SELECT CASE( global%solutFormat )
        CASE ( format_ascii )
          WRITE(fname,'(A,1PE11.5)') &
            trim(global%inDir)//trim(global%casename)//'.plag_stata_', &
            global%timeStamp
          OPEN( if_plag_stats,file=fname,form='formatted',status='old', &
                iostat=errorflag )
      
        CASE ( format_binary )
          WRITE(fname,'(A,1PE11.5)') &
            trim(global%inDir)//trim(global%casename)//'.plag_stat_', &
            global%timeStamp
          OPEN( if_plag_stats,file=fname,form='unformatted',status='old', &
                iostat=errorflag )

        CASE default
          CALL errorstop( global,err_reached_default,__line__ )

      END SELECT ! solutFormat 

    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_file_open,__line__,'File: '//trim(fname) )

    ENDIF ! global%myProcid

! ******************************************************************************
!   Read & broadcast current and integrated time in file, and stats ID 
! ******************************************************************************

    IF ( global%myProcid == masterproc ) THEN
      CALL rflo_readdatafilereal( global,if_plag_stats,global%solutFormat,2,1,rvar )
    ENDIF ! global%myProcid

#ifdef MPI
    CALL mpi_bcast( rvar,2,mpi_rfreal,masterproc,global%mpiComm,global%mpierr )
    IF (global%mpierr /= err_none ) &
      CALL errorstop( global,err_mpi_trouble,__line__ )
#endif

! ==============================================================================
!   Trap error for inconsistent variables in header 
! ==============================================================================
    
    IF ( global%flowType==flow_unsteady .AND. global%currentTime>0._rfreal ) THEN
      IF ( global%currentTime /= rvar(1,1) ) THEN
        WRITE(msg,1000) rvar(1,1),global%currentTime
        CALL errorstop( global,err_time_solution,__line__,msg//' File: '//trim(fname) )
      ENDIF ! currentTime
      
      IF ( global%integrTime  /= rvar(2,1) ) THEN
        WRITE(msg,2000) rvar(2,1),global%integrTime
        CALL errorstop( global,err_time_solution,__line__,msg//' File: '//trim(fname) )
      ENDIF ! integrTime
            
    ENDIF ! global%flowType

! ******************************************************************************
!   Read plagNStat and plagStatId from file
! ******************************************************************************
    
    IF ( global%myProcid == masterproc ) THEN
      IF (global%plagNStat > 0) THEN
        CALL rflo_readdatafileint( global,if_plag_stats,global%solutFormat, &
                                   global%plagNStat+1,1,jvar )
        ntavvar  = jvar(1,1)
        IF ( ntavvar /= global%plagNStat ) THEN
         CALL errorstop( global,err_stats_restart,__line__ )
        END IF ! nTavVar

        plagvarid(1,:) = jvar(2:global%plagNStat+1,1)
        plagvarid(2,:) = mod(plagvarid(1,:),10)
        plagvarid(1,:) = (plagvarid(1,:)-plagvarid(2,:))/10

        DO ind=1,2
        DO l=1,global%plagNStat
          IF ( plagvarid(ind,l) /= global%plagStatId(ind,l) ) &
            CALL errorstop( global,err_stats_restart,__line__ )
        END DO ! l
        END DO ! ind
      ENDIF ! plagNStat

    ENDIF ! myProcid 
      
! ******************************************************************************
!   Read statistics data from all regions 
! ******************************************************************************

    DO ireg=1,global%nRegions

! ==============================================================================
!     Get dimensions and pointers
! ==============================================================================

      ilev = regions(ireg)%currLevel
      CALL rflo_getdimensdummy( regions(ireg),ilev,idcbeg,idcend, &
                                jdcbeg,jdcend,kdcbeg,kdcend )
      CALL rflo_getcelloffset( regions(ireg),ilev,ioff,ijoff )
      ijkbeg = indijk(idcbeg,jdcbeg,kdcbeg,ioff,ijoff)
      ijkend = indijk(idcend,jdcend,kdcend,ioff,ijoff)
      ndim   = ijkend - ijkbeg + 1

      ntav =global%plagNStat

! ==============================================================================
!     Read region number and dimensions (only master)
! ==============================================================================

      IF (global%myProcid == masterproc) THEN
        CALL rflo_readdatafileint( global,if_plag_stats,global%solutFormat,5,1,ivar )
        iregfile  = ivar(1,1)
        ipc       = ivar(2,1)
        jpc       = ivar(3,1)
        kpc       = ivar(4,1)
        ndumcells = ivar(5,1)
      
        IF (iregfile /= ireg) &
          CALL errorstop( global,err_region_number,__line__,'File: '//trim(fname) )
        
        IF ( (ipc /= regions(ireg)%levels(ilev)%grid%ipc) .OR. &
             (jpc /= regions(ireg)%levels(ilev)%grid%jpc) .OR. &
             (kpc /= regions(ireg)%levels(ilev)%grid%kpc)      ) THEN
          WRITE(msg,1005) ireg,ipc,jpc,kpc
          CALL errorstop( global,err_grid_dimensions,__line__,msg )
        ENDIF ! ipc

        IF ( ndumcells /= regions(ireg)%nDumCells ) THEN
          WRITE(msg,1010) ireg,ndumcells,regions(ireg)%nDumCells
          CALL errorstop( global,err_grid_dumcells,__line__,msg )
        ENDIF ! nDumCells

! ==============================================================================
!      Master reads & sends data, others receive them
! ==============================================================================

        ALLOCATE( tavfile(ntav,ndim),stat=errorflag )
        global%error = errorflag
        IF ( global%error /= err_none ) &
          CALL errorstop( global,err_allocate,__line__ )

        CALL rflo_readdatafilereal( global,if_plag_stats,global%solutFormat, &
                                    ntav,ndim,tavfile )

#ifdef MPI
        IF ( regions(ireg)%procid /= masterproc ) THEN
          CALL mpi_send( tavfile,ntav*ndim,mpi_rfreal, &
                         regions(ireg)%procid,ireg,   &
                         global%mpiComm,global%mpierr )
          IF (global%mpierr /= err_none ) &
            CALL errorstop( global,err_mpi_trouble,__line__ )
        ENDIF ! regions(iReg)%procid 
#endif
      ELSE   ! not the master

        IF ( regions(ireg)%procid == global%myProcid ) THEN
          ALLOCATE( tavfile(ntav,ndim),stat=errorflag )
          global%error = errorflag
          IF ( global%error /= err_none ) &
            CALL errorstop( global,err_allocate,__line__ )

#ifdef MPI
          CALL mpi_recv( tavfile,ntav*ndim,mpi_rfreal,masterproc, &
                         ireg,global%mpiComm,status,global%mpierr )
          IF ( global%mpierr /= err_none ) &
            CALL errorstop( global,err_mpi_trouble,__line__ )
#endif
        ENDIF ! regions(iReg)%procid

      ENDIF  ! global%myProcid

! ==============================================================================
!    Copy statistics into data structure
! ==============================================================================
     
      IF ( regions(ireg)%procid == global%myProcid ) THEN
        tav => regions(ireg)%levels(ilev)%plag%tav

        n = 0
        DO k=kdcbeg,kdcend
        DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          n   = n + 1
          ijk = indijk(i,j,k,ioff,ijoff)
          DO l=1,ntav
            tav(l,ijk) = tavfile(l,n)
          ENDDO ! l
        ENDDO ! i
        ENDDO ! j
        ENDDO ! k
      ENDIF ! regions(iReg)%procid

! ==============================================================================
!     Deallocate local array
! ==============================================================================

      IF ( ALLOCATED(tavfile) ) THEN
        DEALLOCATE( tavfile,stat=errorflag )
        global%error = errorflag
        IF ( global%error /= err_none ) &
          CALL errorstop( global,err_deallocate,__line__,'tavFile' )
      ENDIF ! tavFile

    ENDDO     ! iReg

! ******************************************************************************
!   Deallocate temporary data arrays 
! ******************************************************************************

    DEALLOCATE( ivar,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'ivar' )
    
    DEALLOCATE( rvar,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'rvar' )
    
    DEALLOCATE( jvar,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'jvar' )
    
    DEALLOCATE( plagvarid,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'plagVarId' )

! ******************************************************************************
!   Finalize
! ******************************************************************************

    IF ( global%myProcid == masterproc ) THEN
      CLOSE(if_plag_stats,iostat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) &
        CALL errorstop( global,err_file_close,__line__,'File: '//trim(fname) )
    ENDIF

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

! ******************************************************************************
!   Formats 
! ******************************************************************************

1000 FORMAT('Time in file is= ',1pe12.5,' but it should be= ',e12.5,'.')
1005 FORMAT('Region ',i5,', ipc= ',i6,', jpc= ',i6,', kpc= ',i6,'.')
1010 FORMAT('Region ',i5,', # dummy cells=',i2,' but should be= ',i1)
2000 FORMAT('Integration Time in file is= ',1pe12.5,' but it should be= ',e12.5,'.')

  END SUBROUTINE plag_rflo_readstat






!******************************************************************************
!
! Purpose: Wrapper to receive data.
!
! Description: none.
!
! Input: regions = data of all regions
!        iReg    = index of current region.
!
! Output: regions(iReg)%levels%plag%tav = updated plag statitics values
!                                         in dummy cells of current region.
!
! Notes: none.
!
!******************************************************************************

  SUBROUTINE plag_rflo_recvstatbuffwrapper( regions,iReg )

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: regions(:)

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: bctype,ilev,ipatch,ipatchsrc,ireg,iregsrc,npatches 
  
    REAL(KIND=RFREAL), DIMENSION(:,:), POINTER :: tav

    TYPE(t_patch),       POINTER :: ppatch,ppatchsrc
    TYPE(t_region),      POINTER :: pregion,pregionsrc
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(1)%global

    CALL registerfunction( global,'PLAG_RFLO_RecvStatBuffWrapper',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Set dimensions, variables, pointers
! ******************************************************************************    

    pregion => regions(ireg)
    
    ilev     = pregion%currLevel
    npatches = pregion%nPatches

! ******************************************************************************
!   Receive data (regular cells) from other processors 
! ******************************************************************************

    DO ipatch=1,npatches
      ppatch => pregion%levels(ilev)%patches(ipatch)

      bctype    = ppatch%bcType
      iregsrc   = ppatch%srcRegion
      ipatchsrc = ppatch%srcPatch

      IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
          (bctype>=bc_regionint  .AND. bctype<=bc_regionint +bc_range) .OR. &
          (bctype>=bc_regnonconf .AND. bctype<=bc_regnonconf+bc_range) .OR. &
          (bctype>=bc_tra_peri   .AND. bctype<=bc_tra_peri  +bc_range) .OR. &
          (bctype>=bc_rot_peri   .AND. bctype<=bc_rot_peri  +bc_range)) THEN
        IF (regions(iregsrc)%procid /= regions(ireg)%global%myProcid) THEN
          pregionsrc => regions(iregsrc)
          ppatchsrc  => pregionsrc%levels(ilev)%patches(ipatchsrc)

          CALL plag_rflo_recvstatbuff( pregion,pregionsrc,ppatch,ppatchsrc )
        ENDIF ! procid
      
      ENDIF ! bcType 
    ENDDO   ! iPatch

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_recvstatbuffwrapper





!******************************************************************************
!
! Purpose: receive data for dummy cells from adjacent regions being
!          on another processor.
!
! Description: none.
!
! Input: pRegion     = current region
!        pRegionSrc  = source region
!        pPatch      = current patch
!        pPatchSrc   = source patch
!
! Output: pRegion%levels%plag%tav = updated statistics values
!                                   in dummy cells of current region.
!
! Notes: none.
!
!******************************************************************************

  SUBROUTINE plag_rflo_recvstatbuff( pRegion,pRegionSrc,pPatch,pPatchSrc )

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================

    TYPE(t_patch),       POINTER :: ppatch,ppatchsrc
    TYPE(t_region),      POINTER :: pregion,pregionsrc

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: bctype,ilev,ipatch,ipatchsrc,ireg,iregsrc,neqs,npatches
    INTEGER :: ijkbuff,ijkvbuff,ijkd,iv,ivs,lb,n1,n2,ndim
    INTEGER :: i,j,k,idum 
    INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, icoff, ijcoff, source, tag    
    INTEGER :: ivevbeg,ivevend,ivtavbeg,ivtavend,nev,ntav
#ifdef MPI
    INTEGER :: status(mpi_status_size)
#endif

    REAL(KIND=RFREAL), DIMENSION(:,:), POINTER :: ev,tav   

    TYPE(t_buffer_plag), POINTER :: pbuffplag 
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global

    CALL registerfunction( global,'PLAG_RFLO_RecvStatBuff',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Check if the source region is active
! ******************************************************************************
    
    IF ( pregionsrc%active == off ) THEN
      CALL errorstop( global,err_srcregion_off,__line__ )
    ENDIF

! ******************************************************************************
!   Set pointers
! ******************************************************************************

    ilev = pregion%currLevel
    
    ev    => pregion%levels(ilev)%plag%ev
    tav   => pregion%levels(ilev)%plag%tav
    pbuffplag => ppatch%bufferPlag

! ******************************************************************************
!   Set dimensions 
! ******************************************************************************

    CALL rflo_getpatchindices( pregion,ppatch,ilev,ibeg,iend, &
                               jbeg,jend,kbeg,kend )
    CALL rflo_getcelloffset( pregion,ilev,icoff,ijcoff )

    n1   = abs(ppatch%l1end-ppatch%l1beg) + 1   ! here, dimensions of current
    n2   = abs(ppatch%l2end-ppatch%l2beg) + 1   ! and source patch are identical
    ndim = n1*n2*pregion%nDumCells              ! ... but not the # of dummy cells

    ntav = global%plagNStat
    nev  = pregion%levels(ilev)%plag%nEv    
    neqs = nev +ntav       

    ivevbeg = 1
    ivevend = SIZE(ev,dim=1)

    ivtavbeg = 1
    ivtavend = SIZE(tav,dim=1)
 
! ******************************************************************************
!   Receive data
! ******************************************************************************

#ifdef MPI
    source = pregionsrc%procid
    tag    = pregion%localNumber + mpi_patchoff*ppatchsrc%srcPatch
    IF(tag .gt. global%mpiTagMax) tag = mod(tag,global%mpiTagMax)

    CALL mpi_recv( pbuffplag%recvBuffStat,neqs*ndim,mpi_rfreal,   &
                   source,tag,global%mpiComm,status,global%mpierr )
    IF (global%mpierr /= err_none) &
      CALL errorstop( global,err_mpi_trouble,__line__ )
#endif

! ******************************************************************************
!   Copy from buffer to dummy nodes 
! ******************************************************************************
    
    lb       = ppatch%lbound
    ijkbuff  = 0
    ijkvbuff = 0

    DO idum=1,pregion%nDumCells

      SELECT CASE(lb)
      
! ==============================================================================
!       face i=const.
! ==============================================================================

        CASE(1:2)
          IF (lb == 1) i = ibeg - idum
          IF (lb == 2) i = iend + idum
          DO k=kbeg,kend
          DO j=jbeg,jend
            ijkd    = indijk(i,j,k,icoff,ijcoff)
            ijkbuff = ijkbuff + 1
            ivs = 0
            DO iv=ivevbeg,ivevend
              ivs = ivs+1
              ijkvbuff = ijkbuff +(ivs-1)*ndim
              ev(iv,ijkd) = pbuffplag%recvBuffStat(ijkvbuff) 
            ENDDO ! iv
            DO iv=ivtavbeg,ivtavend
              ivs = ivs+1
              ijkvbuff = ijkbuff +(ivs-1)*ndim
              tav(iv,ijkd) = pbuffplag%recvBuffStat(ijkvbuff) 
            ENDDO ! iv
          ENDDO ! j
          ENDDO ! k

! ==============================================================================
!       face j=const.
! ==============================================================================

        CASE(3:4)
          IF (lb == 3) j = jbeg - idum
          IF (lb == 4) j = jend + idum
          DO i=ibeg,iend
          DO k=kbeg,kend
            ijkd    = indijk(i,j,k,icoff,ijcoff)
            ijkbuff = ijkbuff + 1
            ivs = 0
            DO iv=ivevbeg,ivevend
              ivs = ivs+1
              ijkvbuff = ijkbuff +(ivs-1)*ndim
              ev(iv,ijkd) = pbuffplag%recvBuffStat(ijkvbuff) 
            ENDDO ! iv
            DO iv=ivtavbeg,ivtavend
              ivs = ivs+1
              ijkvbuff = ijkbuff +(ivs-1)*ndim
              tav(iv,ijkd) = pbuffplag%recvBuffStat(ijkvbuff) 
            ENDDO ! iv
          ENDDO ! k
          ENDDO ! i

! ==============================================================================
!       face k=const.
! ==============================================================================

        CASE(5:6)
          IF (lb == 5) k = kbeg - idum
          IF (lb == 6) k = kend + idum
          DO j=jbeg,jend
          DO i=ibeg,iend
            ijkd    = indijk(i,j,k,icoff,ijcoff)
            ijkbuff = ijkbuff + 1
            ivs = 0
            DO iv=ivevbeg,ivevend
              ivs = ivs+1
              ijkvbuff = ijkbuff +(ivs-1)*ndim
              ev(iv,ijkd) = pbuffplag%recvBuffStat(ijkvbuff) 
            ENDDO ! iv
            DO iv=ivtavbeg,ivtavend
              ivs = ivs+1
              ijkvbuff = ijkbuff +(ivs-1)*ndim
              tav(iv,ijkd) = pbuffplag%recvBuffStat(ijkvbuff) 
            ENDDO ! iv
          ENDDO  ! i
          ENDDO  ! j
      END SELECT ! lb

    ENDDO ! idum

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_recvstatbuff






! ******************************************************************************
!
! Purpose: Wrapper to send buffer statistics arrays.
!
! Description: None.
!
! Input: regions = dimensions of all regions
!        iReg    = current region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_rflo_sendstatbuffwrapper( regions,iReg )

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: regions(:)

    INTEGER :: ireg

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: bctype,ilev,ipatch,ipatchsrc,iregsrc,npatches

    REAL(RFREAL), POINTER, DIMENSION(:,:) :: tav
  
    TYPE(t_patch),       POINTER :: ppatch
    TYPE(t_region),      POINTER :: pregion,pregionsrc
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(ireg)%global

    CALL registerfunction( global,'PLAG_RFLO_SendStatBuffWrapper',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Load communication buffers
! ******************************************************************************  

    pregion => regions(ireg)
    
    IF ( pregion%procid==global%myProcid .AND. &   ! region active and
         pregion%active==active ) THEN             ! on my processor

! ==============================================================================
!     Set dimensions and pointers
! ==============================================================================

      ilev     = pregion%currLevel
      npatches = pregion%nPatches

! ==============================================================================         
!     Loop over patches 
! ==============================================================================
 
      DO ipatch=1,npatches
        ppatch => pregion%levels(ilev)%patches(ipatch)

        bctype    = ppatch%bcType
        iregsrc   = ppatch%srcRegion
        ipatchsrc = ppatch%srcPatch

        SELECT CASE (bctype)

! ------------------------------------------------------------------------------
!         Conforming region interface
! ------------------------------------------------------------------------------

          CASE( bc_regionconf:bc_regionconf+bc_range )
            pregionsrc => regions(iregsrc)

            IF ( pregionsrc%procid /= global%myProcid ) THEN
              CALL plag_rflo_sendstatbuffconf( pregion,pregionsrc,ppatch )
            ENDIF ! regions(iRegSrc)%procid

! ------------------------------------------------------------------------------
!         Non-conforming region interface (integer)
! ------------------------------------------------------------------------------

          CASE( bc_regionint:bc_regionint+bc_range )
            CALL errorstop( global,err_unknown_bc,__line__ )

! ------------------------------------------------------------------------------
!         Non-conforming region interface (irregular)
! ------------------------------------------------------------------------------
             
          CASE( bc_regnonconf:bc_regnonconf+bc_range)
            CALL errorstop( global,err_unknown_bc,__line__ )           
          
        END select! bcType 

      END DO ! iPatch
    END IF ! pRegion%procid

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_sendstatbuffwrapper






  

! ******************************************************************************
!
! Purpose: Send buffer arrays to dummy cells of the corresponding patch
!          on the adjacent region, residing on another processor.
!
! Description: None.
!
! Input: pRegion    = current region
!        pRegionSrc = source region
!        pPatch     = current patch
!
! Output: None.
!
! Notes: Pertinent for conforming regions.
!
! ******************************************************************************

  SUBROUTINE plag_rflo_sendstatbuffconf( pRegion,pRegionSrc,pPatch )

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
 
    TYPE(t_patch),       POINTER :: ppatch
    TYPE(t_region),      POINTER :: pregion,pregionsrc

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: iv,ivs,irequeststat,npatches,neqs
    INTEGER :: bctype,ipatch,ipatchsrc,iregsrc
    INTEGER :: idum, i, j, k, ijkbuff, ijkvbuff
    INTEGER :: ilev, ibeg, iend, jbeg, jend, kbeg, kend, icoff, ijcoff, ijkc, &
               n1, n2, ndim, dest, tag
    INTEGER :: lb, l1srcdir, l2srcdir, l1beg, l1end, l1step, l2beg, l2end, l2step
    INTEGER :: ivevbeg,ivevend,ivtavbeg,ivtavend,nev,ntav

    LOGICAL :: align
    
    REAL(KIND=RFREAL), DIMENSION(:,:), POINTER :: ev,tav   
    
    TYPE(t_buffer_plag), POINTER :: pbuffplag    
    TYPE(t_plag),        POINTER :: pplag
    TYPE(t_global),      POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global

    CALL registerfunction( global,'PLAG_RFLO_SendStatBuffConf',&
  'PLAG_RFLO_ModStats.F90' )

! ******************************************************************************
!   Check if the source region is active
! ******************************************************************************
    
    IF ( pregionsrc%active == off ) THEN
      CALL errorstop( global,err_srcregion_off,__line__ )
    ENDIF

! ******************************************************************************
!   Set pointers
! ******************************************************************************

    ilev = pregion%currLevel
    
    pplag     => pregion%levels(ilev)%plag
    pbuffplag => ppatch%bufferPlag
    ev        => pplag%ev
    tav       => pplag%tav    

! ******************************************************************************
!   Set dimensions and pointers
! ******************************************************************************
         
    CALL rflo_getpatchindices( pregion,ppatch,ilev,ibeg,iend, &
                               jbeg,jend,kbeg,kend )
    CALL rflo_getcelloffset( pregion,ilev,icoff,ijcoff )

    n1   = abs(ppatch%l1end-ppatch%l1beg) + 1   ! here, dimensions of current
    n2   = abs(ppatch%l2end-ppatch%l2beg) + 1   ! and source patch are identical
    ndim = n1*n2*pregionsrc%nDumCells           ! ... but not the # of dummy cells

    ntav = global%plagNStat
    nev  = pplag%nEv    
    neqs = nev +ntav   
    
    ivevbeg = 1
    ivevend = SIZE(ev,dim=1)
   
    ivtavbeg = 1
    ivtavend = SIZE(tav,dim=1)
 
    irequeststat = pbuffplag%iRequestStat    

! ******************************************************************************
!   Mapping between patches
! ******************************************************************************

                                           l1srcdir =  1
    IF (ppatch%srcL1beg > ppatch%srcL1end) l1srcdir = -1
                                           l2srcdir =  1
    IF (ppatch%srcL2beg > ppatch%srcL2end) l2srcdir = -1

    lb    = ppatch%lbound
    align = ppatch%align

! ******************************************************************************
!   Loop over interior cells of current patch loading buffers 
! ******************************************************************************
    
    ijkbuff  = 0
    ijkvbuff = 0

    DO idum=0,pregionsrc%nDumCells-1

      SELECT CASE(lb)
      
! ==============================================================================
!       face i=const.
! ==============================================================================

        CASE(1:2)
          IF (lb == 1) i = ibeg + idum
          IF (lb == 2) i = iend - idum
          
          IF (align) THEN
            IF (l1srcdir > 0) THEN
              l1beg = jbeg
              l1end = jend
            ELSE
              l1beg = jend
              l1end = jbeg
            ENDIF ! l1SrcDir
            l1step = l1srcdir
            
            IF (l2srcdir > 0) THEN
              l2beg = kbeg
              l2end = kend
            ELSE
              l2beg = kend
              l2end = kbeg
            ENDIF ! l2SrcDir
            l2step = l2srcdir

            DO k=l2beg,l2end,l2step
            DO j=l1beg,l1end,l1step
              ijkc    = indijk(i,j,k,icoff,ijcoff)
              ijkbuff = ijkbuff + 1
              ivs = 0 
              DO iv=ivevbeg,ivevend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = ev(iv,ijkc)
              ENDDO ! iv
              DO iv=ivtavbeg,ivtavend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = tav(iv,ijkc)
              ENDDO ! iv
            ENDDO   ! j
            ENDDO   ! k
          ELSE
            IF (l1srcdir > 0) THEN
              l1beg = kbeg
              l1end = kend
            ELSE
              l1beg = kend
              l1end = kbeg
            ENDIF ! l1SrcDir
            l1step = l1srcdir
        
            IF (l2srcdir > 0) THEN
              l2beg = jbeg
              l2end = jend
            ELSE
              l2beg = jend
              l2end = jbeg
            ENDIF ! l2SrcDir
            l2step = l2srcdir
       
            DO j=l2beg,l2end,l2step
            DO k=l1beg,l1end,l1step
              ijkc    = indijk(i,j,k,icoff,ijcoff)
              ijkbuff = ijkbuff + 1
              ivs = 0
              DO iv=ivevbeg,ivevend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = ev(iv,ijkc)
              ENDDO ! iv
              DO iv=ivtavbeg,ivtavend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = tav(iv,ijkc)
              ENDDO ! iv
            ENDDO   ! j
            ENDDO   ! k
          ENDIF ! align
      
! ==============================================================================
!       face j=const.
! ==============================================================================

        CASE(3:4)
          IF (lb == 3) j = jbeg + idum
          IF (lb == 4) j = jend - idum
          IF (align) THEN
            IF (l1srcdir > 0) THEN
              l1beg = kbeg
              l1end = kend
            ELSE
              l1beg = kend
              l1end = kbeg
            ENDIF ! l1SrcDir
            l1step = l1srcdir
            
            IF (l2srcdir > 0) THEN
              l2beg = ibeg
              l2end = iend
            ELSE
              l2beg = iend
              l2end = ibeg
            ENDIF ! l2SrcDir
            l2step = l2srcdir

            DO i=l2beg,l2end,l2step
            DO k=l1beg,l1end,l1step
              ijkc    = indijk(i,j,k,icoff,ijcoff)
              ijkbuff = ijkbuff + 1
              ivs = 0
              DO iv=ivevbeg,ivevend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = ev(iv,ijkc)
              ENDDO ! iv
              DO iv=ivtavbeg,ivtavend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = tav(iv,ijkc)
              ENDDO ! iv             
            ENDDO ! k
            ENDDO ! i
          ELSE
            IF (l1srcdir > 0) THEN
              l1beg = ibeg
              l1end = iend
            ELSE
              l1beg = iend
              l1end = ibeg
            ENDIF ! l1SrcDir 
            l1step = l1srcdir
            
            IF (l2srcdir > 0) THEN
              l2beg = kbeg
              l2end = kend
            ELSE
              l2beg = kend
              l2end = kbeg
            ENDIF ! l2SrcDir 
            l2step = l2srcdir

            DO k=l2beg,l2end,l2step
            DO i=l1beg,l1end,l1step
              ijkc    = indijk(i,j,k,icoff,ijcoff)
              ijkbuff = ijkbuff + 1
              ivs = 0
              DO iv=ivevbeg,ivevend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = ev(iv,ijkc)
              ENDDO ! iv
              DO iv=ivtavbeg,ivtavend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = tav(iv,ijkc)
              ENDDO ! iv
            ENDDO ! i
            ENDDO ! k
          ENDIF ! align

! ==============================================================================
!       face k=const.
! ==============================================================================

        CASE(5:6)
          IF (lb == 5) k = kbeg + idum
          IF (lb == 6) k = kend - idum
          IF (align) THEN
            IF (l1srcdir > 0) THEN
              l1beg = ibeg
              l1end = iend
            ELSE
              l1beg = iend
              l1end = ibeg
            ENDIF ! l1SrcDir 
            l1step = l1srcdir
            
            IF (l2srcdir > 0) THEN
              l2beg = jbeg
              l2end = jend
            ELSE
              l2beg = jend
              l2end = jbeg
            ENDIF ! l2SrcDir
            l2step = l2srcdir

            DO j=l2beg,l2end,l2step
            DO i=l1beg,l1end,l1step
              ijkc    = indijk(i,j,k,icoff,ijcoff)
              ijkbuff = ijkbuff + 1
              ivs = 0
              DO iv=ivevbeg,ivevend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = ev(iv,ijkc)
              ENDDO ! iv
              DO iv=ivtavbeg,ivtavend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = tav(iv,ijkc)
              ENDDO ! iv
            ENDDO ! i
            ENDDO ! j
          ELSE
            IF (l1srcdir > 0) THEN
              l1beg = jbeg
              l1end = jend
            ELSE
              l1beg = jend
              l1end = jbeg
            ENDIF ! l1SrcDir
            l1step = l1srcdir
            
            IF (l2srcdir > 0) THEN
              l2beg = ibeg
              l2end = iend
            ELSE
              l2beg = iend
              l2end = ibeg
            ENDIF ! l2SrcDir
            l2step = l2srcdir
           
            DO i=l2beg,l2end,l2step
            DO j=l1beg,l1end,l1step
              ijkc    = indijk(i,j,k,icoff,ijcoff)
              ijkbuff = ijkbuff + 1
              ivs = 0
              DO iv=ivevbeg,ivevend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = ev(iv,ijkc)
              ENDDO ! iv
              DO iv=ivtavbeg,ivtavend
                ivs = ivs +1
                ijkvbuff = ijkbuff +(ivs-1)*ndim
                pbuffplag%sendBuffStat(ijkvbuff) = tav(iv,ijkc)
              ENDDO ! iv
            ENDDO ! j
            ENDDO ! i
          ENDIF ! align
      END SELECT ! lb

    ENDDO ! idum

! ******************************************************************************
!   Send data
! ******************************************************************************

#ifdef MPI
    dest = pregionsrc%procid
    tag  = pregionsrc%localNumber + mpi_patchoff*ppatch%srcPatch
    IF(tag .gt. global%mpiTagMax) tag = mod(tag,global%mpiTagMax)

    CALL mpi_isend( pbuffplag%sendBuffStat,neqs*ndim,mpi_rfreal, &
                    dest,tag,global%mpiComm, &
                    pplag%requestsStat(irequeststat),global%mpierr )
  IF (global%mpierr /= err_none) &
    CALL errorstop( global,err_mpi_trouble,__line__ )
#endif

! ******************************************************************************
!   Finalize
! ******************************************************************************

! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_sendstatbuffconf






  



!******************************************************************************
!
! Purpose: Write time averaged statistics of the Lagrangian particles 
!          on Eulerian grid.
!
! Description: the following solution formats are supported:
!              - RocfloMP ASCII
!              - RocfloMP binary
!
! Input: regions = dimensions of all regions
!
! Description: the following solution formats are supported:
!              - RocfloMP ASCII
!              - RocfloMP binary
!
! Input: regions        = dimensions and cons. variables of all regions
!
! Output: to file
!
! Notes: solution is stored only for the current grid level; it is also
!        stored for all dummy cells; all regions are written into one file
!
!******************************************************************************

  SUBROUTINE plag_rflo_writestat(regions)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================

  TYPE(t_region), POINTER :: regions(:)

! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: ireg, i, j, k, l, n

     CHARACTER(2*CHRLEN+17) :: fname

#ifdef MPI
    INTEGER :: status(mpi_status_size)
#endif
    INTEGER :: ilev, iregfile, ipc, jpc, kpc, ndumcells, ndim, ioff, ijoff, ijk
    INTEGER :: idcbeg, jdcbeg, kdcbeg, idcend, jdcend, kdcend, ijkbeg, ijkend
    INTEGER :: errorflag, ntav
    INTEGER, ALLOCATABLE , DIMENSION(:)   :: plagvarid
    INTEGER, ALLOCATABLE , DIMENSION(:,:) :: ivar,jvar

    REAL(RFREAL), POINTER,     DIMENSION(:,:) :: tav
    REAL(RFREAL), ALLOCATABLE, DIMENSION(:,:) :: rvar, tavfile

    TYPE(t_global), POINTER :: global

! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => regions(1)%global

    CALL registerfunction( global,'PLAG_RFLO_WriteStat',&
  'PLAG_RFLO_ModStats.F90' )

    IF (.NOT. global%plagUsed) goto 999

! ******************************************************************************
!   Allocate temporary data arrays 
! ******************************************************************************

    ALLOCATE( ivar(5,1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'ivar' )
    
    ALLOCATE( rvar(2,1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'rvar' )
    
    ALLOCATE( jvar(global%plagNStat+1,1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'rvar' )
    
    ALLOCATE( plagvarid(global%plagNStat+1),stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_allocate,__line__,'plagVarId' ) 

! ******************************************************************************
!   Write verbosity (only master proc.)
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,3X,A)') solver_name,'Writing PLAG statistics solution file...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Open statistics file (only master proc.)
! ******************************************************************************

    IF (global%myProcid == masterproc) THEN

      SELECT CASE( global%solutFormat )
        CASE ( format_ascii )
          WRITE(fname,'(A,1PE11.5)') &
            trim(global%inDir)//trim(global%casename)//'.plag_stata_', &
            global%currentTime
          OPEN( if_plag_stats,file=fname,form='formatted',status='unknown', &
                iostat=errorflag )
      
        CASE ( format_binary )
          WRITE(fname,'(A,1PE11.5)') &
            trim(global%inDir)//trim(global%casename)//'.plag_stat_', &
            global%currentTime
          OPEN( if_plag_stats,file=fname,form='unformatted',status='unknown', &
                iostat=errorflag )
        
        CASE default
          CALL errorstop( global,err_reached_default,__line__ )

      END SELECT ! solutFormat 

      global%error = errorflag
      IF ( global%error /= err_none ) &
        CALL errorstop( global,err_file_open,__line__,'File: '//trim(fname) )

! ******************************************************************************
!     Write current and integrated time to file 
! ******************************************************************************

      rvar(1,1) = global%currentTime
      rvar(2,1) = global%integrTime

      CALL rflo_writedatafilereal( global,if_plag_stats,global%solutFormat,&
                                   2,1,rvar )

! ******************************************************************************
!     Write plagNStat and plagStatId to file
! ******************************************************************************
      
      IF ( global%plagNStat > 0 ) THEN
        jvar(1,1)    = global%plagNStat
        plagvarid(:) = global%plagStatId(1,:)*10 + global%plagStatId(2,:)
        jvar(2:global%plagNStat+1,1) = plagvarid(1:global%plagNStat)
        CALL rflo_writedatafileint( global,if_plag_stats,global%solutFormat, &
                                    global%plagNStat+1,1,jvar )
      ENDIF ! plagNStat

    ENDIF ! global%myProcid

! ******************************************************************************
!   Write statistics data 
! ******************************************************************************

    DO ireg=1,global%nRegions

! ==============================================================================
!     Get dimensions and pointers
! ==============================================================================

      ilev = regions(ireg)%currLevel
      CALL rflo_getdimensdummy( regions(ireg),ilev,idcbeg,idcend, &
                                jdcbeg,jdcend,kdcbeg,kdcend )
      CALL rflo_getcelloffset( regions(ireg),ilev,ioff,ijoff )
      ijkbeg = indijk(idcbeg,jdcbeg,kdcbeg,ioff,ijoff)
      ijkend = indijk(idcend,jdcend,kdcend,ioff,ijoff)
      ndim   = ijkend - ijkbeg + 1

      ntav = global%plagNStat

! ==============================================================================
!    Allocate memory for data field
! ==============================================================================

      IF ( regions(ireg)%procid==global%myProcid .OR. &
           global%myProcid==masterproc                ) THEN
        ALLOCATE( tavfile(ntav,ndim),stat=errorflag )
        global%error = errorflag
        IF ( global%error /= err_none ) &
          CALL errorstop( global,err_allocate,__line__,'tavFile' )
      ENDIF ! regions(iReg)%procid

! ==============================================================================
!     Copy statistics into data structure
! ==============================================================================

      IF ( regions(ireg)%procid == global%myProcid ) THEN
        tav => regions(ireg)%levels(ilev)%plag%tav

        n = 0
        DO k=kdcbeg,kdcend
        DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          n   = n + 1
          ijk = indijk(i,j,k,ioff,ijoff)
          DO l=1,ntav
            tavfile(l,n) = tav(l,ijk)
          ENDDO ! l
        ENDDO ! i
        ENDDO ! j
        ENDDO ! k
      ENDIF ! regions(iReg)%procid 

! ==============================================================================
!     Write region number and dimensions (only master)
! ==============================================================================

      IF ( global%myProcid == masterproc ) THEN
        ivar(1,1) = ireg
        ivar(2,1) = regions(ireg)%levels(ilev)%grid%ipc
        ivar(3,1) = regions(ireg)%levels(ilev)%grid%jpc
        ivar(4,1) = regions(ireg)%levels(ilev)%grid%kpc
        ivar(5,1) = regions(ireg)%nDumCells
        CALL rflo_writedatafileint( global,if_plag_stats,global%solutFormat,5,1,ivar )
      ENDIF ! global%myProcid

! ==============================================================================
!     Master receives and writes data, others send them
! ==============================================================================

      IF (global%myProcid == masterproc) THEN
#ifdef MPI
        IF ( regions(ireg)%procid /= masterproc ) THEN
          CALL mpi_recv( tavfile,ntav*ndim,mpi_rfreal, &
                         regions(ireg)%procid,ireg, &
                         global%mpiComm,status,global%mpierr )
          IF ( global%mpierr /= err_none ) &
            CALL errorstop( global,err_mpi_trouble,__line__ )
        ENDIF ! regions(iReg)%procid
#endif
        CALL rflo_writedatafilereal( global,if_plag_stats,global%solutFormat, &
                                     ntav,ndim,tavfile )

      ELSE   ! not the master
#ifdef MPI
        IF (regions(ireg)%procid == global%myProcid) THEN
          CALL mpi_send( tavfile,ntav*ndim,mpi_rfreal,masterproc, &
                         ireg,global%mpiComm,global%mpierr )
          IF ( global%mpierr /= err_none ) &
            CALL errorstop( global,err_mpi_trouble,__line__ )
      ENDIF
#endif
      ENDIF  ! global%myProcid

! ==============================================================================
!     Deallocate local array
! ==============================================================================

      IF (ALLOCATED(tavfile)) THEN
        DEALLOCATE( tavfile,stat=errorflag )
        global%error = errorflag
        IF ( global%error /= err_none ) &
          CALL errorstop( global,err_deallocate,__line__,'tavFile' )
      ENDIF

    ENDDO     ! iReg

! ******************************************************************************
!   Deallocate temporary data arrays 
! ******************************************************************************

    DEALLOCATE( ivar,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'ivar' )
    
    DEALLOCATE( rvar,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'rvar' )
   
    DEALLOCATE( jvar,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'jvar' )
   
    DEALLOCATE( plagvarid,stat=errorflag )
    global%error = errorflag
    IF ( global%error /= err_none ) &
      CALL errorstop( global,err_deallocate,__line__,'plagVarId' )

! ******************************************************************************
!   Write verbosity (only master proc.)
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN
      WRITE(stdout,'(A,3X,A)') solver_name,'Writing PLAG statistics solution file done...'
    END IF ! global%verbLevel

! ******************************************************************************
!   Finalize
! ******************************************************************************

    IF ( global%myProcid == masterproc ) THEN
      CLOSE(if_plag_stats,iostat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) &
        CALL errorstop( global,err_file_close,__line__,'File: '//trim(fname) )
    ENDIF

! ******************************************************************************
!   End  
! ******************************************************************************

999  CONTINUE
    CALL deregisterfunction( global )

  END SUBROUTINE plag_rflo_writestat







END MODULE plag_rflo_modstats

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: PLAG_RFLO_ModStats.F90,v $
! Revision 1.9  2009/03/02 00:19:36  mtcampbe
! Added some ifdefs around Rocflo to disable particle injection on INFLOW
! boundaries and added some checks around MPI tags utilizing a new global
! data item, global%mpiTagMax.
!
! Revision 1.8  2008/12/06 08:44:34  mtcampbe
! Updated license.
!
! Revision 1.7  2008/11/19 22:17:47  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.6  2006/04/07 15:19:24  haselbac
! Removed tabs
!
! Revision 1.5  2005/06/19 07:13:57  wasistho
! refrained from communicating and writing statistics when plag is not used
!
! Revision 1.4  2005/03/08 00:59:00  fnajjar
! Bug fix for ijkvBuff with multiple assignments
!
! Revision 1.3  2005/03/07 17:37:37  fnajjar
! Bug fix in loading iRequestStat
!
! Revision 1.2  2005/02/16 14:47:41  fnajjar
! Added infrastructure for on-processor copying and MPI-based communication
!
! Revision 1.1  2005/01/08 20:44:32  fnajjar
! Initial import for PLAG statistics
!
! ******************************************************************************