PLAG_ModSurfStats.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 surface statistics.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: PLAG_ModSurfStats.F90,v 1.13 2008/12/06 08:44:34 mtcampbe Exp $
!
! Copyright: (c) 2004 by the University of Illinois
!
! ******************************************************************************

MODULE plag_modsurfstats

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

  USE modbuildfilenames, ONLY: buildfilenamesteady, &
                               buildfilenameunsteady

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

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

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

! ==============================================================================
! Public functions
! ==============================================================================
  
  PUBLIC :: plag_createsurfstats, & 
            plag_destroysurfstats, & 
            plag_gathersurfstats, &
            plag_readsurfstatswrapper, &
            plag_writesurfstatswrapper

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

  PRIVATE :: plag_closesurfstats, & 
             plag_createsurfstatskernel, &
             plag_destroysurfstatskernel, &
             plag_initsurfstats, &
             plag_nullifysurfstats, &
             plag_opensurfstatsascii, &
             plag_opensurfstatsbinary, &
             plag_readsurfstatsascii, & 
             plag_readsurfstatsbinary, & 
             plag_readsurfstatskernelascii, & 
             plag_readsurfstatskernelbinary, & 
             plag_setbinindex, &
             plag_writesurfstatsascii, & 
             plag_writesurfstatsbinary, &
             plag_writesurfstatskernelascii, & 
             plag_writesurfstatskernelbinary

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

  CONTAINS
  
  
  
  
  
! *******************************************************************************
!
! Purpose: Close surface statitstics file.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_closesurfstats(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    INTEGER :: errorflag,ifile
    TYPE(t_global), POINTER :: global    
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_CloseSurfStats',&
  'PLAG_ModSurfStats.F90')
             
    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Closing PLAG surface statistics file...'
    END IF ! global%verbLevel             

    ifile = if_plag_surf_stats
        
! ******************************************************************************
!   Close file
! ******************************************************************************

    CLOSE(ifile,iostat=errorflag)
    global%error = errorflag      
    IF ( global%error /= err_none ) THEN 
      CALL errorstop(global,err_file_close,__line__)
    END IF ! global%error

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_closesurfstats  
  
  
  
  
  

! *******************************************************************************
!
! Purpose: Create surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_createsurfstats(pRegion)

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

    INTEGER :: errorflag,ipatch
    TYPE(t_global), POINTER :: global
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************
    
    global => pregion%global
    
    CALL registerfunction(global,'PLAG_CreateSurfStats',&
  'PLAG_ModSurfStats.F90')
    
    pgrid => pregion%grid
        
! ******************************************************************************
!   Allocate memory
! ******************************************************************************
    
! ==============================================================================
!   Patch data
! ==============================================================================

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN        
        ALLOCATE(ppatch%statsPlag(ppatch%nBFaces),stat=errorflag)
        global%error = errorflag
        IF ( global%error /= err_none ) THEN 
          CALL errorstop(global,err_allocate,__line__,'pPatch%statsPlag')
        END IF ! global%error
     
        CALL plag_createsurfstatskernel(pregion,ppatch%statsPlag)                           
      END IF ! pPatch%plotStatsFlag
    END DO ! iPatch
    
! ==============================================================================
!   Plane data
! ==============================================================================

! ******************************************************************************
!   Initialize memory
! ******************************************************************************
    
! ==============================================================================
!   Patch data
! ==============================================================================

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)
      IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
        CALL plag_initsurfstats(pregion,ppatch%statsPlag)                           
      END IF ! pPatch%plotStatsFlag
    END DO ! iPatch
    
! ==============================================================================
!   Plane data
! ==============================================================================

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_createsurfstats
  
  
  
  
  

! *******************************************************************************
!
! Purpose: Kernel to create surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_createsurfstatskernel(pRegion,pStatsPlag)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: pregion
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
 
! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: errorflag,ifl,nbfaces,nbins,ncont,nvars

    TYPE(t_global), POINTER :: global    

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

    global => pregion%global
 
    CALL registerfunction(global,'PLAG_CreateSurfStatsKernel',&
  'PLAG_ModSurfStats.F90')

    nbfaces = SIZE(pstatsplag,dim=1)   

    ncont = pregion%plagInput%nCont
    nvars = plag_surf_stats_last +ncont
     
! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY

! ******************************************************************************
!   Allocate memory
! ******************************************************************************

    DO ifl = 1,nbfaces
      ALLOCATE(pstatsplag(ifl)%nHits(nbins),stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_allocate,__line__,'pStatsPlag%nHits')
      END IF ! global%error
        
      ALLOCATE(pstatsplag(ifl)%vars(nvars,nbins),stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_allocate,__line__,'pStatsPlag%vars')
      END IF ! global%error     
    END DO ! ifl                           

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_createsurfstatskernel 






! *******************************************************************************
!
! Purpose: Destroy surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_destroysurfstats(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    INTEGER :: errorflag,ipatch
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_DestroySurfStats',&
  'PLAG_ModSurfStats.F90')
    
    pgrid => pregion%grid
        
! ******************************************************************************
!   Deallocate memory
! ******************************************************************************
        
    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN       
        CALL plag_destroysurfstatskernel(global,ppatch%statsPlag)

        DEALLOCATE(ppatch%statsPlag,stat=errorflag)
        global%error = errorflag
        IF ( global%error /= err_none ) THEN 
          CALL errorstop(global,err_deallocate,__line__,'pPatch%statsPlag')
        END IF ! global%error
      END IF ! pPatch%plotStatsFlag  
    END DO ! iPatch

! ******************************************************************************
!   Nullify memory
! ******************************************************************************
    
    CALL plag_nullifysurfstats(pregion)    

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_destroysurfstats
  





! *******************************************************************************
!
! Purpose: Kernel to destroy surface statistics.
!
! Description: None.
!
! Input:
!   global              Pointer to global
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_destroysurfstatskernel(global,pStatsPlag)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_global), POINTER :: global  
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
 
! ==============================================================================
!   Locals
! ==============================================================================
    
    INTEGER :: errorflag,ifl,nbfaces

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

    CALL registerfunction(global,'PLAG_DestroySurfStatsKernel',&
  'PLAG_ModSurfStats.F90')

    nbfaces = SIZE(pstatsplag,dim=1)        

! ******************************************************************************
!   Deallocate memory
! ******************************************************************************
        
    DO ifl = 1,nbfaces
      DEALLOCATE(pstatsplag(ifl)%nHits,stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_deallocate,__line__,'pStatsPlag%nHits')
      END IF ! global%error
        
      DEALLOCATE(pstatsplag(ifl)%vars,stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_deallocate,__line__,'pPatch%statsPlag%vars')
      END IF ! global%error     
    END DO ! ifl                           

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_destroysurfstatskernel
  
  
  
  
  

! *******************************************************************************
!
! Purpose: Gather surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pPlag               Pointer to plag
!   pStatsPlag          Pointer to statsPlag
!   ifl                 Local face index  
!   iPcl                Particle index
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_gathersurfstats(pRegion,pPlag,pStatsPlag,ifl,iPcl,thetaAngle)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion   
    TYPE(t_plag), POINTER :: pplag
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag

    INTEGER,INTENT(IN) :: ifl,ipcl
    REAL(KIND=RFREAL),INTENT(IN) :: thetaangle
    
! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: binmethod,errorflag,ibin,icont,nbins,ncont
    REAL(KIND=RFREAL) :: diambin,diamminmicron,diammaxmicron,diamplagmicron, &
                         enerplag,massplag,momemagnplag

    TYPE(t_global), POINTER :: global     

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

    global => pregion%global

    CALL registerfunction(global,'PLAG_GatherSurfStats',&
  'PLAG_ModSurfStats.F90')

    ncont = pregion%plagInput%nCont

! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY

! TEMPORARY
  binmethod = bin_method_linear
! binMethod = pRegion%plagInput%binMethod
! END TEMPORARY

! TEMPORARY
  diamminmicron =  10.0_rfreal
  diammaxmicron = 210.0_rfreal
  diambin       = (diammaxmicron-diamminmicron)/REAL(nbins,kind=rfreal)
! END TEMPORARY
  
  diamplagmicron = pplag%dv(dv_plag_diam,ipcl)*1.0e+6_rfreal       
  
  massplag = sum( pplag%cv(pplag%cvPlagMass(:),ipcl) )
  
  momemagnplag = sqrt( pplag%cv(cv_plag_xmom,ipcl)**2.0_rfreal &
                      +pplag%cv(cv_plag_ymom,ipcl)**2.0_rfreal &
                      +pplag%cv(cv_plag_zmom,ipcl)**2.0_rfreal )
  
  enerplag = pplag%cv(cv_plag_ener,ipcl)

! ******************************************************************************
!   Gather statistics
! ******************************************************************************
    
! ==============================================================================
!   Determine bin
! ==============================================================================

    ibin = plag_setbinindex(global,binmethod,nbins,diamplagmicron,diambin, &
                            diamminmicron,diammaxmicron)

! ==============================================================================
!   Cummulate data on bins
! ==============================================================================

    pstatsplag(ifl)%nHits(ibin) = pstatsplag(ifl)%nHits(ibin) +1
 
    pstatsplag(ifl)%vars(plag_surf_stats_diam3,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_diam3,ibin) +diamplagmicron**3.0_rfreal 

    pstatsplag(ifl)%vars(plag_surf_stats_diam4,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_diam4,ibin) +diamplagmicron**4.0_rfreal 

    pstatsplag(ifl)%vars(plag_surf_stats_theta,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_theta,ibin) +thetaangle 

    pstatsplag(ifl)%vars(plag_surf_stats_mome1,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_mome1,ibin) +momemagnplag

    pstatsplag(ifl)%vars(plag_surf_stats_mome2,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_mome2,ibin) +momemagnplag**2.0_rfreal

    pstatsplag(ifl)%vars(plag_surf_stats_ener,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_ener,ibin) +enerplag

    pstatsplag(ifl)%vars(plag_surf_stats_mass,ibin) = &
      pstatsplag(ifl)%vars(plag_surf_stats_mass,ibin) +massplag

    DO icont = 1, ncont
      pstatsplag(ifl)%vars(plag_surf_stats_last+icont,ibin) = &
       pstatsplag(ifl)%vars(plag_surf_stats_last+icont,ibin)  &
      +pplag%cv(pplag%cvPlagMass(icont),ipcl)
    END DO ! iCont 

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_gathersurfstats
  





! *******************************************************************************
!
! Purpose: Initialize surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_initsurfstats(pRegion,pStatsPlag)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
     
    TYPE(t_region), POINTER :: pregion  
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
  
! ==============================================================================
!   Locals
! ==============================================================================
 
    INTEGER :: errorflag,ibin,ifl,ipatch,ivar,nbfaces,nbins,ncont,nvars

    TYPE(t_global), POINTER :: global  

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

    global => pregion%global
 
    CALL registerfunction(global,'PLAG_InitSurfStats',&
  'PLAG_ModSurfStats.F90')

    nbfaces = SIZE(pstatsplag,dim=1)          

    ncont = pregion%plagInput%nCont
    nvars = plag_surf_stats_last +ncont

! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY 
       
! ******************************************************************************
!   Initialize memory
! ******************************************************************************
        
    DO ifl = 1,nbfaces
      DO ibin = 1, nbins
        pstatsplag(ifl)%nHits(ibin) = 0
      END DO ! iBin

      DO ibin = 1, nbins
      DO ivar = 1, nvars
        pstatsplag(ifl)%vars(ivar,ibin) = 0.0_rfreal
      END DO ! iVar
      END DO ! iBin                  
    END DO ! ifl                      

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_initsurfstats






! *******************************************************************************
!
! Purpose: Nullify surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_nullifysurfstats(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    INTEGER :: errorflag,ipatch
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_DestroySurfStats',&
  'PLAG_ModSurfStats.F90')
    
    pgrid => pregion%grid

! ******************************************************************************
!   Nullify pointer
! ******************************************************************************
     
! ==============================================================================
!     Patch data
! ==============================================================================

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
        nullify(ppatch%statsPlag)
      END IF ! pPatch%plotStatsFlag
    END DO ! iPatch
   
! ==============================================================================
!     Plane data
! ==============================================================================

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_nullifysurfstats






! *******************************************************************************
!
! Purpose: Open surface statistics file in ASCII format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   fileStatus          File status
!
! Output: 
!   fileExists          Flag indicating existence of file
!
! Notes:
!   1. File may not exist when it is opened in postprocessor before solver 
!      was run, so need to deal with this gracefully.
!
! ******************************************************************************

  SUBROUTINE plag_opensurfstatsascii(pRegion,fileStatus,fileExists)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    INTEGER, INTENT(IN) :: filestatus
    LOGICAL, INTENT(OUT), OPTIONAL :: fileexists
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename
    INTEGER :: errorflag,ifile
    TYPE(t_global), POINTER :: global 
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_OpenSurfStatsASCII',&
  'PLAG_ModSurfStats.F90')
         
    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Opening ASCII PLAG surface statistics file...'
    END IF ! global%verbLevel     
     
    ifile = if_plag_surf_stats
   
! ******************************************************************************
!   Build file name
! ******************************************************************************

    CALL buildfilenameunsteady(global,filedest_outdir,'.plag_ssta', & 
                               pregion%iRegionGlobal,global%currentTime, & 
                               ifilename)

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN                                             
      WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', & 
                                       pregion%iRegionGlobal
      WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name,'Current time:', & 
                                          global%currentTime  
    END IF ! global%verbLevel

! ******************************************************************************
!   Open file
! ******************************************************************************

    IF ( filestatus == file_status_old ) THEN 
      INQUIRE(file=ifilename,exist=fileexists)
    
      IF ( fileexists .EQV. .true. ) THEN     
        OPEN(ifile,file=ifilename,form="FORMATTED",status="OLD", &
             iostat=errorflag)
        global%error = errorflag          
        IF ( global%error /= err_none ) THEN 
          CALL errorstop(global,err_file_open,__line__,ifilename)
        END IF ! global%error
      END IF ! fileExists
    ELSE IF ( filestatus == file_status_unknown ) THEN 
      OPEN(ifile,file=ifilename,form="FORMATTED",status="UNKNOWN", &
           iostat=errorflag)
      global%error = errorflag          
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_file_open,__line__,ifilename)
      END IF ! global%error      
    ELSE 
      CALL errorstop(global,err_reached_default,__line__)      
    END IF ! fileStatus       

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_opensurfstatsascii






! *******************************************************************************
!
! Purpose: Open surface statistics file in binary format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   fileStatus          File status
!
! Output: 
!   fileExists          Flag indicating existence of file
!
! Notes: 
!   1. File may not exist when it is opened in postprocessor before solver 
!      was run, so need to deal with this gracefully.
!
! ******************************************************************************

  SUBROUTINE plag_opensurfstatsbinary(pRegion,fileStatus,fileExists)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    INTEGER, INTENT(IN) :: filestatus
    LOGICAL, INTENT(OUT), OPTIONAL :: fileexists
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename
    INTEGER :: errorflag,ifile
    TYPE(t_global), POINTER :: global 
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_OpenSurfStatsBinary',&
  'PLAG_ModSurfStats.F90')
         
    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Opening binary PLAG surface statistics file...'
    END IF ! global%verbLevel     
     
    ifile = if_plag_surf_stats     
     
! ******************************************************************************
!   Build file name
! ******************************************************************************

    CALL buildfilenameunsteady(global,filedest_outdir,'.plag_sst', & 
                               pregion%iRegionGlobal,global%currentTime, & 
                               ifilename)

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN                                             
      WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', & 
                                       pregion%iRegionGlobal
      WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name,'Current time:', & 
                                          global%currentTime  
    END IF ! global%verbLevel

! ******************************************************************************
!   Open file
! ******************************************************************************

    IF ( filestatus == file_status_old ) THEN 
      INQUIRE(file=ifilename,exist=fileexists)
    
      IF ( fileexists .EQV. .true. ) THEN     
        OPEN(ifile,file=ifilename,form="UNFORMATTED",status="OLD", &
             iostat=errorflag)
        global%error = errorflag          
        IF ( global%error /= err_none ) THEN 
          CALL errorstop(global,err_file_open,__line__,ifilename)
        END IF ! global%error
      END IF ! fileExists
    ELSE IF ( filestatus == file_status_unknown ) THEN 
      OPEN(ifile,file=ifilename,form="UNFORMATTED",status="UNKNOWN", &
           iostat=errorflag)
      global%error = errorflag          
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_file_open,__line__,ifilename)
      END IF ! global%error      
    ELSE 
      CALL errorstop(global,err_reached_default,__line__)      
    END IF ! fileStatus         

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_opensurfstatsbinary
  
  
  
  
  
  
! *******************************************************************************
!
! Purpose: Read surface statistics in ASCII format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_readsurfstatsascii(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename,sectionstring
    INTEGER :: errorflag,ifile,ifl,ipatch,loopcounter
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_ReadSurfStatsASCII',&
  'PLAG_ModSurfStats.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Reading ASCII PLAG surface statistics file...'
    END IF ! global%verbLevel
    
    pgrid => pregion%grid

    ifile = if_plag_surf_stats
         
! ******************************************************************************
!   Header and general information
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_low ) THEN  
      WRITE(stdout,'(A,3X,A)') solver_name,'Header information...'
    END IF ! global%verbLevel

    READ(ifile,'(A)') sectionstring
    IF ( trim(sectionstring) /= '# ROCFLU surface statistics file' ) THEN 
      CALL errorstop(global,err_invalid_marker,__line__,sectionstring) 
    END IF ! TRIM

! ******************************************************************************
!   Read rest of file
! ******************************************************************************

    loopcounter = 0

    DO ! set up infinite loop
      loopcounter = loopcounter + 1

      READ(ifile,'(A)') sectionstring

      SELECT CASE ( trim(sectionstring) ) 

! ==============================================================================
!       Patch data
! ==============================================================================            

        CASE ( '# Patch' )
          IF ( global%myProcid == masterproc .AND. &
               global%verbLevel > verbose_low ) THEN  
            WRITE(stdout,'(A,3X,A)') solver_name,'Patch...'
          END IF ! global%verbLevel 

         DO ipatch = 1,pgrid%nPatches
           ppatch => pregion%patches(ipatch)
           IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
             CALL plag_readsurfstatskernelascii(pregion,ppatch%statsPlag)
           END IF ! pPatch%plotStatsFlag
         END DO ! iPatch

! ==============================================================================
!       Plane data
! ==============================================================================            

! ==============================================================================
!       End marker
! ==============================================================================
      
        CASE ( '# End' ) 
          IF ( global%myProcid == masterproc .AND. & 
               global%verbLevel > verbose_low ) THEN  
            WRITE(stdout,'(A,3X,A)') solver_name,'End marker...'
          END IF ! global%verbLevel           

          EXIT
      
! ==============================================================================
!       Invalid section string
! ==============================================================================
      
        CASE default
          IF ( global%verbLevel > verbose_low ) THEN  
            WRITE(stdout,'(A,3X,A)') solver_name,sectionstring
          END IF ! verbosityLevel           

          CALL errorstop(global,err_invalid_marker,__line__,sectionstring)                               
      END SELECT ! TRIM
  
! ==============================================================================
!     Guard against infinite loop - unnecessary because of read errors?
! ==============================================================================
  
      IF ( loopcounter >= limit_infinite_loop ) THEN 
        CALL errorstop(global,err_infinite_loop,__line__)
      END IF ! loopCounter
    END DO ! <empty>       

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_readsurfstatsascii






! *******************************************************************************
!
! Purpose: Read surface statistics in binary format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_readsurfstatsbinary(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename,sectionstring
    INTEGER :: errorflag,ifile,ifl,ipatch,loopcounter
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_ReadSurfStatsBinary',&
  'PLAG_ModSurfStats.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Reading ASCII PLAG surface statistics file...'
    END IF ! global%verbLevel
    
    pgrid => pregion%grid

    ifile = if_plag_surf_stats
         
! ******************************************************************************
!   Header and general information
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_low ) THEN  
      WRITE(stdout,'(A,3X,A)') solver_name,'Header information...'
    END IF ! global%verbLevel

    READ(ifile) sectionstring
    IF ( trim(sectionstring) /= '# ROCFLU surface statistics file' ) THEN 
      CALL errorstop(global,err_invalid_marker,__line__,sectionstring) 
    END IF ! TRIM

! ******************************************************************************
!   Read rest of file
! ******************************************************************************

    loopcounter = 0

    DO ! set up infinite loop
      loopcounter = loopcounter + 1

      READ(ifile) sectionstring

      SELECT CASE ( trim(sectionstring) ) 

! ==============================================================================
!       Patch data
! ==============================================================================            

        CASE ( '# Patch' )
          IF ( global%myProcid == masterproc .AND. &
               global%verbLevel > verbose_low ) THEN  
            WRITE(stdout,'(A,3X,A)') solver_name,'Patch...'
          END IF ! global%verbLevel 

         DO ipatch = 1,pgrid%nPatches
           ppatch => pregion%patches(ipatch)
           IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
             CALL plag_readsurfstatskernelbinary(pregion,ppatch%statsPlag)
           END IF ! pPatch%plotStatsFlag
         END DO ! iPatch

! ==============================================================================
!       Plane data
! ==============================================================================            

! ==============================================================================
!       End marker
! ==============================================================================
      
        CASE ( '# End' ) 
          IF ( global%myProcid == masterproc .AND. & 
               global%verbLevel > verbose_low ) THEN  
            WRITE(stdout,'(A,3X,A)') solver_name,'End marker...'
          END IF ! global%verbLevel           

          EXIT
      
! ==============================================================================
!       Invalid section string
! ==============================================================================
      
        CASE default
          IF ( global%verbLevel > verbose_low ) THEN  
            WRITE(stdout,'(A,3X,A)') solver_name,sectionstring
          END IF ! verbosityLevel           

          CALL errorstop(global,err_invalid_marker,__line__,sectionstring)                               
      END SELECT ! TRIM
  
! ==============================================================================
!     Guard against infinite loop - unnecessary because of read errors?
! ==============================================================================
  
      IF ( loopcounter >= limit_infinite_loop ) THEN 
        CALL errorstop(global,err_infinite_loop,__line__)
      END IF ! loopCounter
    END DO ! <empty>       

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_readsurfstatsbinary






! *******************************************************************************
!
! Purpose: Kernel to read surface statistics in ASCII format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_readsurfstatskernelascii(pRegion,pStatsPlag)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
   
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename,sectionstring
    INTEGER :: errorflag,ibin,ifile,ifl,ipatch,ivar,nbfaces,nbins,ncont,nvars  

    TYPE(t_global), POINTER :: global     

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

    global => pregion%global

    CALL registerfunction(global,'PLAG_ReadSurfStatsKernelASCII',&
  'PLAG_ModSurfStats.F90')

    ifile = if_plag_surf_stats

    nbfaces = SIZE(pstatsplag,dim=1)

    ncont = pregion%plagInput%nCont
    nvars = plag_surf_stats_last +ncont
        
! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY

! ******************************************************************************
!   Read surface statistics to file
! ******************************************************************************
    
! ==============================================================================
!   Number of hits
! ==============================================================================

    READ(ifile,'(8(I10))') ((pstatsplag(ifl)%nHits(ibin),ibin=1,nbins),&
                             ifl=1,nbfaces)        

! ==============================================================================
!   Surface variables
! ==============================================================================

    READ(ifile,'(5(E23.16))') (((pstatsplag(ifl)%vars(ivar,ibin),         &
                                 ivar=1,nvars),ibin=1,nbins),ifl=1,nbfaces)
   
! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction(global)

  END SUBROUTINE plag_readsurfstatskernelascii






! *******************************************************************************
!
! Purpose: Kernel to read surface statistics in binary format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_readsurfstatskernelbinary(pRegion,pStatsPlag)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
   
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename,sectionstring
    INTEGER :: errorflag,ibin,ifile,ifl,ipatch,ivar,nbfaces,nbins,ncont,nvars  

    TYPE(t_global), POINTER :: global     

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

    global => pregion%global

    CALL registerfunction(global,'PLAG_ReadSurfStatsKernelBinary',&
  'PLAG_ModSurfStats.F90')

    ifile = if_plag_surf_stats

    nbfaces = SIZE(pstatsplag,dim=1)

    ncont = pregion%plagInput%nCont
    nvars = plag_surf_stats_last +ncont
        
! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY

! ******************************************************************************
!   Read surface statistics to file
! ******************************************************************************
    
! ==============================================================================
!   Number of hits
! ==============================================================================

    READ(ifile) ((pstatsplag(ifl)%nHits(ibin),ibin=1,nbins),ifl=1,nbfaces)        

! ==============================================================================
!   Surface variables
! ==============================================================================

    READ(ifile) (((pstatsplag(ifl)%vars(ivar,ibin),         &
                   ivar=1,nvars),ibin=1,nbins),ifl=1,nbfaces)
   
! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction(global)

  END SUBROUTINE plag_readsurfstatskernelbinary












! *******************************************************************************
!
! Purpose: Wrapper for reading surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_readsurfstatswrapper(pRegion)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================

    LOGICAL :: fileexists    
    TYPE(t_global), POINTER :: global  
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_ReadSurfStatsWrapper',&
  'PLAG_ModSurfStats.F90')
    
    SELECT CASE(global%solutFormat)
      CASE( format_ascii )
        CALL plag_opensurfstatsascii(pregion,file_status_old,fileexists)
      
      IF ( fileexists .EQV. .true. ) THEN
        CALL plag_readsurfstatsascii(pregion)
        CALL plag_closesurfstats(pregion) 
      END IF ! fileExists

      CASE( format_binary )
        CALL plag_opensurfstatsbinary(pregion,file_status_old,fileexists)
      
         IF ( fileexists .EQV. .true. ) THEN 
           CALL plag_readsurfstatsbinary(pregion)
           CALL plag_closesurfstats(pregion)      
         END IF ! fileExists
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! global%solutFormat   

! ******************************************************************************
!   Write warning if file is missing
! ******************************************************************************
    
    IF ( fileexists .EQV. .false. ) THEN 
      global%warnCounter = global%warnCounter + 1
    
      IF ( global%myProcid == masterproc .AND. & 
           global%verbLevel > verbose_none ) THEN 
        WRITE(stdout,'(A,2(1X,A))') solver_name,'*** WARNING ***', &
                              'Plag surface statistics file missing, not read.'        
      END IF ! global%myProcid
    END IF ! fileExists     

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_readsurfstatswrapper
  
  
  
  
  

! *******************************************************************************
!
! Purpose: Set bin index based on method selected..
!
! Description: None.
!
! Input:
!   global              Pointer to global
!   binMethod           Binning method selected
!   diamMicron          Particle diameter in microns
!   diamMin             Minimum diameter
!   diamMax             Maximum diameter
!
! Output: Bin index.
!
! Notes: None.
!
! ******************************************************************************

  INTEGER FUNCTION plag_setbinindex(global,binMethod,nBins,diamMicron,&
                                    diambin,diammin,diammax)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_global), POINTER :: global   
    
    INTEGER, INTENT(IN) :: binmethod,nbins
    REAL(KIND=RFREAL), INTENT(IN) :: diambin,diammicron,diammax,diammin
    
! ==============================================================================
!   Locals
! ==============================================================================

    INTEGER :: errorflag

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

    CALL registerfunction(global,'PLAG_SetBinMethod',&
  'PLAG_ModSurfStats.F90')

! ******************************************************************************
!   Determine bin index
! ******************************************************************************

    SELECT CASE(binmethod)
      CASE( bin_method_linear)
        plag_setbinindex = min(int((diammicron-diammin)/diambin)+1,nbins)

      CASE default 
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! binMethod

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

    CALL deregisterfunction(global)

  END FUNCTION plag_setbinindex  
  
  
  
  

  
! *******************************************************************************
!
! Purpose: Write surface statistics in ASCII format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_writesurfstatsascii(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename,sectionstring
    INTEGER :: errorflag,ifile,ifl,ipatch
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_WriteSurfStatsASCII',&
  'PLAG_ModSurfStats.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Writing ASCII PLAG surface statistics file...'
    END IF ! global%verbLevel
    
    pgrid => pregion%grid

    ifile = if_plag_surf_stats
         
! ******************************************************************************
!   Header and general information
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_low ) THEN  
      WRITE(stdout,'(A,3X,A)') solver_name,'Header information...'
    END IF ! global%verbLevel

    sectionstring = '# ROCFLU surface statistics file'
    WRITE(ifile,'(A)') trim(sectionstring)
                            
! ******************************************************************************
!   Write patch data
! ******************************************************************************

    sectionstring = '# Patch'
    WRITE(ifile,'(A)') trim(sectionstring)
        
    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
        CALL plag_writesurfstatskernelascii(pregion,ppatch%statsPlag)
      END IF ! pPatch%plotStatsFlag
    END DO ! iPatch
    
! ******************************************************************************
!   Write plane data
! ******************************************************************************

! ******************************************************************************
!   End marker
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_low ) THEN  
      WRITE(stdout,'(A,3X,A)') solver_name,'End marker...'
    END IF ! global%verbLevel

    sectionstring = '# End'
    WRITE(ifile,'(A)') trim(sectionstring)

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_writesurfstatsascii





! *******************************************************************************
!
! Purpose: Write surface statistics in binary format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_writesurfstatsbinary(pRegion)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: sectionstring
    INTEGER :: errorflag,ifile,ipatch
    TYPE(t_global), POINTER :: global    
    TYPE(t_grid), POINTER :: pgrid
    TYPE(t_patch), POINTER :: ppatch
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_WriteSurfStatsBinary',&
  'PLAG_ModSurfStats.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_none ) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name, & 
                               'Writing binary PLAG surface statistics file...'
    END IF ! global%verbLevel
    
    pgrid => pregion%grid

    ifile = if_plag_surf_stats

! ******************************************************************************
!   Header and general information
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_low ) THEN  
      WRITE(stdout,'(A,3X,A)') solver_name,'Header information...'
    END IF ! global%verbLevel

    sectionstring = '# ROCFLU surface statistics file'
    WRITE(ifile) sectionstring
              
! ******************************************************************************
!   Write patch data
! ******************************************************************************

    sectionstring = '# Patch'
    WRITE(ifile) sectionstring
        
    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( ppatch%plotStatsFlag .EQV. .true. ) THEN
        CALL plag_writesurfstatskernelbinary(pregion,ppatch%statsPlag)
      END IF ! pPatch%plotStatsFlag
    END DO ! iPatch
    
! ******************************************************************************
!   Write plane data
! ******************************************************************************
    
          

! ******************************************************************************
!   End marker
! ******************************************************************************

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel > verbose_low ) THEN  
      WRITE(stdout,'(A,3X,A)') solver_name,'End marker...'
    END IF ! global%verbLevel

    sectionstring = '# End'
    WRITE(ifile) sectionstring

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

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_writesurfstatsbinary






! *******************************************************************************
!
! Purpose: Kernel to write surface statistics in ASCII format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_writesurfstatskernelascii(pRegion,pStatsPlag)

    IMPLICIT NONE
        
! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
   
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename,sectionstring
    INTEGER :: errorflag,ibin,ifile,ifl,ipatch,ivar,nbfaces,nbins,ncont,nvars  

    TYPE(t_global), POINTER :: global     

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

    global => pregion%global

    CALL registerfunction(global,'PLAG_WriteSurfStatsKernelASCII',&
  'PLAG_ModSurfStats.F90')

    ifile = if_plag_surf_stats

    nbfaces = SIZE(pstatsplag,dim=1)

    ncont = pregion%plagInput%nCont
    nvars = plag_surf_stats_last +ncont
      
! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY

! ******************************************************************************
!   Write surface statistics to file
! ******************************************************************************
    
! ==============================================================================
!   Number of hits
! ==============================================================================

    WRITE(ifile,'(8(I10))') ((pstatsplag(ifl)%nHits(ibin),ibin=1,nbins),&
                              ifl=1,nbfaces)        

! ==============================================================================
!   Surface variables
! ==============================================================================

    WRITE(ifile,'(5(E23.16))') (((pstatsplag(ifl)%vars(ivar,ibin),         &
                                  ivar=1,nvars),ibin=1,nbins),ifl=1,nbfaces)
   
! ******************************************************************************
!   End  
! ******************************************************************************

    CALL deregisterfunction(global)

  END SUBROUTINE plag_writesurfstatskernelascii






! *******************************************************************************
!
! Purpose: Kernel to write surface statistics in binary format.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!   pStatsPlag          Pointer to statsPlag
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_writesurfstatskernelbinary(pRegion,pStatsPlag)

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

    TYPE(t_region), POINTER :: pregion   

    TYPE(t_surfstats_plag), DIMENSION(:), POINTER :: pstatsplag
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    CHARACTER(CHRLEN) :: ifilename
    INTEGER :: errorflag,ibin,ifile,ifl,ipatch,ivar,nbfaces,nbins,ncont,nvars

    TYPE(t_global), POINTER :: global  
 
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global

    CALL registerfunction(global,'PLAG_WriteSurfStatsKernelBinary',&
  'PLAG_ModSurfStats.F90')

    ifile = if_plag_surf_stats 

    nbfaces = SIZE(pstatsplag,dim=1)

    ncont = pregion%plagInput%nCont
    nvars = plag_surf_stats_last +ncont
          
! TEMPORARY
    nbins = 20
!    nBins = pRegion%plagInput%nBins
! END TEMPORARY

! ******************************************************************************
!   Write surface statistics to file
! ******************************************************************************
    
! ==============================================================================
!   Number of hits
! ==============================================================================

    WRITE(ifile) ((pstatsplag(ifl)%nHits(ibin),ibin=1,nbins),ifl=1,nbfaces)       
    
! ==============================================================================
!   Surface variables
! ==============================================================================

    WRITE(ifile) (((pstatsplag(ifl)%vars(ivar,ibin),         &
                    ivar=1,nvars),ibin=1,nbins),ifl=1,nbfaces)

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

    CALL deregisterfunction(global)

  END SUBROUTINE plag_writesurfstatskernelbinary





! *******************************************************************************
!
! Purpose: Wrapper for writing surface statistics.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE plag_writesurfstatswrapper(pRegion)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************
    
! ==============================================================================
!   Arguments
! ==============================================================================
    
    TYPE(t_region), POINTER :: pregion
    
! ==============================================================================
!   Locals
! ==============================================================================
    
    TYPE(t_global), POINTER :: global  
    
! ******************************************************************************
!   Start, set pointers and variables
! ******************************************************************************

    global => pregion%global
    
    CALL registerfunction(global,'PLAG_WriteSurfStatsWrapper',&
  'PLAG_ModSurfStats.F90')
    
    SELECT CASE( global%solutFormat ) 
      CASE( format_ascii )
        CALL plag_opensurfstatsascii(pregion,file_status_unknown)
        CALL plag_writesurfstatsascii(pregion)
        CALL plag_closesurfstats(pregion)            
      CASE( format_binary )
        CALL plag_opensurfstatsbinary(pregion,file_status_unknown)
        CALL plag_writesurfstatsbinary(pregion)
        CALL plag_closesurfstats(pregion)      
      CASE default 
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! global%solutFormat   
    
! ******************************************************************************
!   End
! ******************************************************************************

    CALL deregisterfunction(global)

  END SUBROUTINE plag_writesurfstatswrapper






END MODULE plag_modsurfstats

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: PLAG_ModSurfStats.F90,v $
! Revision 1.13  2008/12/06 08:44:34  mtcampbe
! Updated license.
!
! Revision 1.12  2008/11/19 22:17:47  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.11  2006/05/12 22:50:32  fnajjar
! Fixed binary write of strings by remove TRIM
!
! Revision 1.10  2006/05/11 18:17:44  fnajjar
! Fixed inconsistent naming for register functions
!
! Revision 1.9  2006/05/09 14:38:10  fnajjar
! Bug for a formatted read from a binary file
!
! Revision 1.8  2006/05/04 13:54:58  fnajjar
! Added and activated binary read for patch statistics file
!
! Revision 1.7  2006/05/02 17:46:23  fnajjar
! Allowed surface statistics to be gathered on active patches
!
! Revision 1.6  2006/04/07 15:19:23  haselbac
! Removed tabs
!
! Revision 1.5  2005/12/30 16:27:18  fnajjar
! Changed diamMaxMicron to have size of 10 in bins
!
! Revision 1.4  2005/01/06 22:26:43  fnajjar
! Implemented clean version for PLAG_SetBinIndex
!
! Revision 1.3  2005/01/06 16:00:31  fnajjar
! Reformulated IDINT to INT with KIND in PLAG_SetBinIndex
!
! Revision 1.2  2005/01/05 19:26:53  fnajjar
! Fixed definition of diamBin to maintain bounds
!
! Revision 1.1  2004/12/21 15:07:02  fnajjar
! Initial import of surface statistics module
!
! ******************************************************************************