source: palm/trunk/UTIL/read_prt_data.f90 @ 1771

 PROGRAM read_prt_data
!--------------------------------------------------------------------------------!
! This file is part of PALM.
!
! PALM is free software: you can redistribute it and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2014  Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
! Initial revision.
!
! Former revisions:
! -----------------
! $Id: analyze_particle_data.f90 1310 2014-03-14 08:01:56Z raasch $
!
! Description:
! ------------
! This routine reads the particle data generated by PALM, and enables user
! analysis. Compile and execute this program in the folder where your particle 
! data (_000000, _000001, ...) is stored.
!-------------------------------------------------------------------------------!

    IMPLICIT NONE

    CHARACTER(LEN=7)    ::  i_proc_char
    CHARACTER (LEN=80)  ::  rtext
    CHARACTER (LEN=110) ::  run_description_header

    REAL(KIND=8) ::  simulated_time
    INTEGER      ::  ip, i_proc=0, i_proc_end, jp, kp,                         &
                     max_number_of_particle_groups, number_of_particles,       &
                     number_of_particle_groups, n, nxl,                        &
                     nxr, nys, nyn, nzb, nzt, nbgp, status

    INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::  prt_count

    TYPE particle_type
       SEQUENCE
       REAL(KIND=8) ::  radius, age, age_m, dt_sum, dvrp_psize, e_m,           &
                        origin_x, origin_y, origin_z, rvar1, rvar2, rvar3,     &
                        speed_x, speed_y, speed_z, weight_factor, x, y, z
      INTEGER       ::  class, group, tailpoints, tail_id
      LOGICAL       ::  particle_mask
      INTEGER       ::  block_nr
    END TYPE particle_type

    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  particles_temp
    TYPE(particle_type), DIMENSION(:), POINTER     ::  particles

    TYPE  grid_particle_def
       TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  particles
    END TYPE grid_particle_def

    TYPE(grid_particle_def), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  grid_particles

    TYPE particle_groups_type
        SEQUENCE
        REAL(KIND=8) ::  density_ratio, radius, exp_arg, exp_term
    END TYPE particle_groups_type

    TYPE(particle_groups_type), DIMENSION(:), ALLOCATABLE ::  particle_groups

    LOGICAL ::  found
!
!-- Check if file from PE0 exists and terminate program if it doesn't.
    WRITE (i_proc_char,'(''_'',I6.6)')  i_proc
    INQUIRE ( FILE=i_proc_char, EXIST=found )
!
!-- Estimate the number of files (equal to the number of PEs which
!-- have been used in PALM)
    DO  WHILE ( found )
       OPEN ( 85, FILE=i_proc_char, FORM='UNFORMATTED' )
       i_proc = i_proc + 1
       WRITE (i_proc_char,'(''_'',I6.6)')  i_proc
       INQUIRE ( FILE=i_proc_char, EXIST=found )
       CLOSE( 85 )
    ENDDO
!
!-- Info-output
    PRINT*, ''
    PRINT*, '*** read_prt ***'
    IF ( i_proc /= 0 )  THEN
       PRINT*, '***', i_proc, ' file(s) found'
    ELSE
       PRINT*, '+++ file _000000 not found'
       PRINT*, '    program terminated'
       STOP
    ENDIF
!
!-- Set number of files and opens them sequentially
    i_proc_end = i_proc-1

    DO  i_proc = 0, i_proc_end
!
!--    Open particle data file for each process
       WRITE (i_proc_char,'(''_'',I6.6)')  i_proc
       OPEN ( 85, FILE=i_proc_char, FORM='UNFORMATTED' ) 
!
!--    Read general description of file and inform user
       READ ( 85 )  run_description_header
       READ ( 85 )  rtext
       IF ( i_proc == 0 )  THEN
          PRINT*, ' '
          PRINT*, '*** data description header:'
          PRINT*, '*** ', run_description_header
          PRINT*, '*** ', rtext
          PRINT*, ' '
       ENDIF
       PRINT*, '*** data of file', i_proc, 'is analysed'
       READ ( 85 )  number_of_particle_groups, max_number_of_particle_groups

       IF ( .NOT. ALLOCATED(particle_groups) )  THEN
          ALLOCATE( particle_groups(max_number_of_particle_groups) )
       ENDIF
    
       READ ( 85 )  particle_groups
       READ ( 85 )  nxl, nxr, nys, nyn, nzb, nzt, nbgp

       IF ( ALLOCATED(prt_count)  .OR.  ALLOCATED(grid_particles) )  THEN
          DEALLOCATE( prt_count, grid_particles )
       ENDIF

       ALLOCATE( prt_count(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp),     &
                 grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )

!
!--    Start individual time loop
       DO
          READ( 85, IOSTAT=status)  simulated_time
          IF ( status < 0 )  THEN
             PRINT*, '*** end of file reached'
             EXIT
          ENDIF
          PRINT*, '*** time of analyzed data set:', simulated_time

          READ ( 85 )  prt_count

          DO  ip = nxl, nxr
             DO  jp = nys, nyn
                DO  kp = nzb+1, nzt
                   number_of_particles = prt_count(kp,jp,ip)
                   IF ( number_of_particles <= 0 )  CYCLE
                   ALLOCATE( grid_particles(kp,jp,ip)%particles(1:number_of_particles) )
                   ALLOCATE( particles_temp(1:number_of_particles) )
                   READ ( 85 )  particles_temp
                   grid_particles(kp,jp,ip)%particles = particles_temp
                   DEALLOCATE( particles_temp )
                ENDDO
             ENDDO
          ENDDO
!
!--       This part can be used for user analysis
          DO  ip = nxl, nxr
             DO  jp = nys, nyn
                DO  kp = nzb+1, nzt
                   number_of_particles = prt_count(kp,jp,ip)
                   IF ( number_of_particles <= 0 )  CYCLE
                   particles => grid_particles(kp,jp,ip)%particles
!
!--                Add your analysis here
!                   DO  n = 1, number_of_particles
!
!                   ENDDO

                ENDDO
             ENDDO
          ENDDO
!
!--       Deallocate grid_particles%particles in case of more than one timestep
          DO  ip = nxl, nxr
             DO  jp = nys, nyn
                DO  kp = nzb+1, nzt
                   number_of_particles = prt_count(kp,jp,ip)
                   IF ( number_of_particles <= 0 )  CYCLE
                   DEALLOCATE( grid_particles(kp,jp,ip)%particles )
                ENDDO
             ENDDO
          ENDDO     

       ENDDO ! end of time loop

       CLOSE( 85 )

    ENDDO ! end of file loop

 END PROGRAM read_prt_data