PROGRAM read_prt_data !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! 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 . ! ! Copyright 1997-2018 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: read_prt_data.f90 2696 2017-12-14 17:12:51Z kanani $ ! Adapted for changed particle attributes ! ! 2696 2017-12-14 17:12:51Z kanani ! Corrected "Former revisions" section ! ! 2696 2017-12-14 17:12:51Z kanani ! Change in file header (GPL part) ! ! 2123 2017-01-18 12:49:59Z hoffmann ! Particle variables renamed: tailpoints, tail_id, and dvrp_psize to id1, id2, ! and user ! ! 1771 2016-02-29 10:57:56Z hoffmann ! Initial revision. ! ! 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) :: aux1 !< auxiliary multi-purpose feature REAL(KIND=8) :: aux2 !< auxiliary multi-purpose feature REAL(KIND=8) :: radius !< radius of particle REAL(KIND=8) :: age !< age of particle REAL(KIND=8) :: age_m !< REAL(KIND=8) :: dt_sum !< REAL(KIND=8) :: e_m !< interpolated sgs tke REAL(KIND=8) :: origin_x !< origin x-position of particle (changed cyclic bc) REAL(KIND=8) :: origin_y !< origin y-position of particle (changed cyclic bc) REAL(KIND=8) :: origin_z !< origin z-position of particle (changed cyclic bc) REAL(KIND=8) :: rvar1 !< REAL(KIND=8) :: rvar2 !< REAL(KIND=8) :: rvar3 !< REAL(KIND=8) :: speed_x !< speed of particle in x REAL(KIND=8) :: speed_y !< speed of particle in y REAL(KIND=8) :: speed_z !< speed of particle in z REAL(KIND=8) :: weight_factor !< weighting factor REAL(KIND=8) :: x !< x-position REAL(KIND=8) :: y !< y-position REAL(KIND=8) :: z !< z-position INTEGER(KIND=4) :: class !< radius class needed for collision INTEGER(KIND=4) :: group !< number of particle group INTEGER(KIND=8) :: id !< particle ID (64 bit integer) LOGICAL :: particle_mask !< if this parameter is set to false the particle will be deleted INTEGER(KIND=4) :: block_nr !< number for sorting (removable?) 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