Index: palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4731)
+++ palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4778)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! output of particle time series added
+!
+! 4731 2020-10-07 13:25:11Z schwenkel
 ! Move exchange_horiz from time_integration to modules
 !
@@ -209,4 +212,12 @@
     USE cpulog,                                                                                    &
         ONLY:  cpu_log, log_point, log_point_s
+
+    USE data_output_particle_mod,                                                                  &
+        ONLY:  dop_active,                                                                         &
+               dop_finalize,                                                                       &
+               dop_init,                                                                           &
+               dop_output_tseries,                                                                 &
+               dop_prt_axis_dimension,                                                             &
+               dop_last_active_particle
 
     USE indices,                                                                                   &
@@ -266,8 +277,4 @@
                surf_lsm_h,                                                                         &
                surf_usm_h
-
-!-- Next lines are in preparation for the output of particle data
-!    USE data_output_particle_mod,                                              &
-!        ONLY:  dop_active, dop_init, dop_output_tseries
 
 #if defined( __parallel )
@@ -326,4 +333,5 @@
     INTEGER(iwp) ::  trnp_count_sum                               !< parameter for particle exchange of PEs
 
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  id_counter !< number of particles initialized in each grid box
     INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array_particles   !< sequence of random array for particle
 
@@ -410,20 +418,21 @@
     PRIVATE
 
-    PUBLIC lpm_parin,                                                                              &
-           lpm_header,                                                                             &
-           lpm_init_arrays,                                                                        &
-           lpm_init,                                                                               &
-           lpm_actions,                                                                            &
+    PUBLIC lagrangian_particle_model
+
+    PUBLIC lpm_actions,                                                                            &
+           lpm_check_parameters,                                                                   &
            lpm_data_output_ptseries,                                                               &
            lpm_exchange_horiz_bounds,                                                              &
+           lpm_header,                                                                             &
+           lpm_init,                                                                               &
            lpm_interaction_droplets_ptq,                                                           &
+           lpm_init_arrays,                                                                        &
+           lpm_last_actions,                                                                       &
+           lpm_parin,                                                                              &
+           lpm_rrd_global,                                                                         &
+           lpm_rrd_local,                                                                          &
            lpm_rrd_local_particles,                                                                &
-           lpm_wrd_local,                                                                          &
-           lpm_rrd_global,                                                                         &
            lpm_wrd_global,                                                                         &
-           lpm_rrd_local,                                                                          &
-           lpm_check_parameters
-
-    PUBLIC lagrangian_particle_model
+           lpm_wrd_local
 
     INTERFACE lpm_check_parameters
@@ -533,4 +542,8 @@
        MODULE PROCEDURE dealloc_particles_array
     END INTERFACE dealloc_particles_array
+
+    INTERFACE lpm_last_actions
+       MODULE PROCEDURE lpm_last_actions
+    END INTERFACE lpm_last_actions
 
     INTERFACE lpm_sort_and_delete
@@ -629,4 +642,5 @@
        collision_kernel,                                                                           &
        curvature_solution_effects,                                                                 &
+       data_output_pts,                                                                            &
        deallocate_memory,                                                                          &
        density_ratio,                                                                              &
@@ -651,7 +665,7 @@
        particle_advection_interpolation,                                                           &
        particle_maximum_age,                                                                       &
-       part_output,                                                                                &
-       part_inc,                                                                                   &
-       part_percent,                                                                               &
+       pts_id_file,                                                                                &
+       pts_increment,                                                                              &
+       pts_percentage,                                                                             &
        pdx,                                                                                        &
        pdy,                                                                                        &
@@ -861,4 +875,6 @@
  SUBROUTINE lpm_check_parameters
 
+    LOGICAL ::  id_file_exists = .FALSE.
+
 !
 !-- Collision kernels:
@@ -880,4 +896,5 @@
 
     END SELECT
+
     IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
 
@@ -897,7 +914,33 @@
        message_string = 'nested runs in combination with cloud droplets ' //                       &
                         'is not implemented'
-          CALL message( 'lpm_check_parameters', 'PA0687', 1, 2, 0, 6, 0 )
+       CALL message( 'lpm_check_parameters', 'PA0687', 1, 2, 0, 6, 0 )
     ENDIF
 
+    IF ( pts_increment > 1  .AND.  pts_percentage < 100.0_wp )  THEN
+       message_string = 'pts_increment and pts_percentage cannot be set both '
+       CALL message( 'lpm_check_parameters', 'PA0...', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( pts_increment < 1 )  THEN
+       message_string = 'pts_increment must be > 1'
+       CALL message( 'lpm_check_parameters', 'PA0...', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( pts_percentage > 100.0_wp )  THEN
+       message_string = 'pts_percentage must be < 100'
+       CALL message( 'lpm_check_parameters', 'PA0...', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    INQUIRE( FILE = pts_id_file, EXIST = id_file_exists )
+#if defined( __netcdf4_parallel )
+!
+!-- Check if individual particles timeseries is set
+    IF ( pts_increment  > 1  .AND.  dt_dopts /= 9999999.9_wp  .OR.                                 &
+         pts_percentage < 100.0_wp  .AND.  dt_dopts /= 9999999.9_wp  .OR.                          &
+         id_file_exists  .AND.  dt_dopts /= 9999999.9_wp  )                                        &
+    THEN
+       dop_active = .TRUE.
+    ENDIF
+#endif
 
  END SUBROUTINE lpm_check_parameters
@@ -948,4 +991,6 @@
     INTEGER(iwp) ::  j                           !<
     INTEGER(iwp) ::  k                           !<
+
+    LOGICAL  ::  read_restart                    !<
 
     REAL(wp) ::  div                             !<
@@ -1229,13 +1274,17 @@
 !-- For the first model run of a possible job chain initialize the particles, otherwise read the
 !-- particle data from restart file.
+    read_restart = .FALSE.
     IF ( TRIM( initializing_actions ) == 'read_restart_data'                                       &
          .AND.  read_particles_from_restartfile )  THEN
        CALL lpm_rrd_local_particles
+
+       read_restart = .TRUE.
     ELSE
 !
 !--    Allocate particle arrays and set attributes of the initial set of particles, which can be
 !--    also periodically released at later times.
-       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                         &
-                 grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
+       ALLOCATE( grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr),                                        &
+                 id_counter(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                        &
+                 prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
        number_of_particles = 0
@@ -1243,5 +1292,5 @@
 !
 !--    Initialize counter for particle IDs
-       grid_particles%id_counter = 1
+       id_counter = 1
 !
 !--    Initialize all particles with dummy values (otherwise errors may occur within restart runs).
@@ -1314,6 +1363,8 @@
 
 !
-!-- next line is in preparation for particle data output
-!    CALL dop_init
+!-- Output of particle time series
+    IF ( dop_active )  THEN
+       CALL dop_init( read_restart )
+    ENDIF
 
 !
@@ -1520,5 +1571,5 @@
     prt_count   = local_count
 !
-!-- Calculate particle IDs
+!-- Calculate particle IDs (for new particles only)
     DO  ip = nxl, nxr
        DO  jp = nys, nyn
@@ -1527,15 +1578,13 @@
              IF ( number_of_particles <= 0 )  CYCLE
              particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
-
-             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
-
-                particles(n)%id = 10000_idp**3 * grid_particles(kp,jp,ip)%id_counter +             &
-                                  10000_idp**2 * kp + 10000_idp * jp + ip
+             DO  n = local_start(kp,jp,ip), number_of_particles
+
+                particles(n)%id = 10000_idp**3 * id_counter(kp,jp,ip) + 10000_idp**2 * kp +        &
+                                  10000_idp * jp + ip
 !
 !--             Count the number of particles that have been released before
-                grid_particles(kp,jp,ip)%id_counter = grid_particles(kp,jp,ip)%id_counter + 1
+                id_counter(kp,jp,ip) = id_counter(kp,jp,ip) + 1
 
              ENDDO
-
           ENDDO
        ENDDO
@@ -2682,4 +2731,8 @@
     CALL cpu_log( log_point(36), 'data_output_ptseries', 'start' )
 
+    IF ( dop_active )  THEN
+       CALL dop_output_tseries
+    ENDIF
+
     IF ( myid == 0 )  THEN
 !
@@ -3052,8 +3105,9 @@
        FLUSH(9)
 
-       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                         &
-                 grid_particles(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
-       ALLOCATE( prt_global_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE( grid_particles(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                    &
+                 id_counter(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                        &
+                 prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                         &
+                 prt_global_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
 !
 !--    Open restart file for read, if not already open, and do not allow usage of shared-memory I/O
@@ -3068,4 +3122,5 @@
 
        prt_count = 0
+       CALL rrd_mpi_io( 'id_counter', id_counter )
        CALL rrd_mpi_io( 'prt_count', prt_count )
        CALL rrd_mpi_io( 'prt_global_index', prt_global_index )
@@ -3147,10 +3202,18 @@
     LOGICAL, INTENT(OUT)  ::  found
 
-    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::  tmp_3d   !<
-
+    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::  tmp_3d          !<
+    INTEGER(iwp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::  tmp_3d_int  !<
 
     found = .TRUE.
 
     SELECT CASE ( restart_string(1:length) )
+
+        CASE ( 'id_counter' )
+           IF ( .NOT. ALLOCATED( id_counter ) )  THEN
+              ALLOCATE( id_counter(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+           ENDIF
+           IF ( k == 1 )  READ ( 13 )  tmp_3d_int
+           id_counter(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                 &
+              tmp_3d_int(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
         CASE ( 'pc_av' )
@@ -3349,4 +3412,8 @@
 #endif
 
+       IF ( ALLOCATED( id_counter ) )  THEN
+          CALL wrd_write_string( 'id_counter' )
+          WRITE ( 14 )  id_counter
+       ENDIF
 
        IF ( ALLOCATED( seq_random_array_particles ) )  THEN
@@ -3355,7 +3422,5 @@
        ENDIF
 
-
     ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
-
 
        IF ( ALLOCATED( seq_random_array_particles ) )  THEN
@@ -3396,5 +3461,6 @@
        ENDIF
 
-       CALL wrd_mpi_io( 'prt_count', prt_count )
+       CALL wrd_mpi_io( 'id_counter', id_counter )
+       CALL wrd_mpi_io( 'prt_count',  prt_count )
 
        start_index = start_index_on_pe
@@ -3437,4 +3503,10 @@
        WRITE ( 14 )  curvature_solution_effects
 
+       CALL wrd_write_string( 'dop_last_active_particle' )
+       WRITE ( 14 )  dop_last_active_particle
+
+       CALL wrd_write_string( 'dop_prt_axis_dimension' )
+       WRITE ( 14 )  dop_prt_axis_dimension
+
        CALL wrd_write_string( 'interpolation_simple_corrector' )
        WRITE ( 14 )  interpolation_simple_corrector
@@ -3459,4 +3531,6 @@
        CALL wrd_mpi_io( 'bc_par_ns', bc_par_ns )
        CALL wrd_mpi_io( 'bc_par_t', bc_par_t )
+       CALL wrd_mpi_io( 'dop_last_active_particle', dop_last_active_particle )
+       CALL wrd_mpi_io( 'dop_prt_axis_dimension', dop_prt_axis_dimension )
        CALL wrd_mpi_io( 'last_particle_release_time', last_particle_release_time )
        CALL wrd_mpi_io( 'number_of_particle_groups', number_of_particle_groups )
@@ -3508,4 +3582,10 @@
           READ ( 13 )  curvature_solution_effects
 
+       CASE ( 'dop_last_active_particle' )
+          READ ( 13 )  dop_last_active_particle
+
+       CASE ( 'dop_prt_axis_dimension' )
+          READ ( 13 )  dop_prt_axis_dimension
+
        CASE ( 'interpolation_simple_corrector' )
           READ ( 13 )  interpolation_simple_corrector
@@ -3551,4 +3631,6 @@
     CALL rrd_mpi_io( 'bc_par_ns', bc_par_ns )
     CALL rrd_mpi_io( 'bc_par_t', bc_par_t )
+    CALL rrd_mpi_io( 'dop_prt_axis_dimension', dop_prt_axis_dimension )
+    CALL rrd_mpi_io( 'dop_last_active_particle', dop_last_active_particle )
     CALL rrd_mpi_io( 'last_particle_release_time', last_particle_release_time )
     CALL rrd_mpi_io( 'number_of_particle_groups', number_of_particle_groups )
@@ -3577,4 +3659,20 @@
 
  END SUBROUTINE lpm_rrd_global_mpi
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Last actions before PALM finishes.
+!------------------------------------------------------------------------------!
+ SUBROUTINE lpm_last_actions
+
+!
+!-- Close NETCDF file for individual particle timeseries
+    IF ( dop_active )  THEN
+       CALL dop_finalize
+    ENDIF
+
+ END SUBROUTINE lpm_last_actions
 
 
