Changeset 1359


Ignore:
Timestamp:
Apr 11, 2014 5:15:14 PM (8 years ago)
Author:
hoffmann
Message:

new Lagrangian particle structure integrated

Location:
palm/trunk/SOURCE
Files:
1 added
2 deleted
38 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1338 r1359  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# mod_particle_attributes added, lpm_sort_arrays removed,
     23# lpm_extend_particle_array removed
    2324#
    2425# Former revisions:
     
    174175        lpm_data_output_particles.f90 lpm_droplet_collision.f90 \
    175176        lpm_droplet_condensation.f90 lpm_exchange_horiz.f90 \
    176         lpm_extend_particle_array.f90 lpm_extend_tails.f90 \
     177        lpm_extend_tails.f90 \
    177178        lpm_extend_tail_array.f90 lpm_init.f90 lpm_init_sgs_tke.f90 \
    178179        lpm_pack_arrays.f90 lpm_read_restart_file.f90 lpm_release_set.f90 \
    179         lpm_set_attributes.f90 lpm_sort_arrays.f90 \
     180        lpm_set_attributes.f90 \
    180181        lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \
    181182        ls_forcing.f90 message.f90 microphysics.f90 modules.f90 mod_kinds.f90 \
    182         netcdf.f90 nudging.f90 \
     183        mod_particle_attributes.f90 netcdf.f90 nudging.f90 \
    183184        package_parin.f90 palm.f90 parin.f90 plant_canopy_model.f90 poisfft.f90 \
    184185        poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
     
    261262data_log.o: modules.o mod_kinds.o
    262263data_output_dvrp.o: modules.o cpulog.o mod_kinds.o
    263 data_output_mask.o: modules.o cpulog.o mod_kinds.o
     264data_output_mask.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
    264265data_output_profiles.o: modules.o cpulog.o mod_kinds.o
    265 data_output_ptseries.o: modules.o cpulog.o mod_kinds.o
     266data_output_ptseries.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
    266267data_output_spectra.o: modules.o cpulog.o mod_kinds.o
    267268data_output_tseries.o: modules.o cpulog.o mod_kinds.o
    268 data_output_2d.o: modules.o cpulog.o mod_kinds.o
    269 data_output_3d.o: modules.o cpulog.o mod_kinds.o
     269data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
     270data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
    270271diffusion_e.o: modules.o mod_kinds.o
    271272diffusion_s.o: modules.o mod_kinds.o
     
    286287inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
    287288init_1d_model.o: modules.o mod_kinds.o
    288 init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o ls_forcing.o
     289init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o \
     290        ls_forcing.o lpm_init.o
    289291init_advec.o: modules.o mod_kinds.o
    290292init_cloud_physics.o: modules.o mod_kinds.o
     
    304306local_tremain.o: modules.o cpulog.o mod_kinds.o
    305307local_tremain_ini.o: modules.o cpulog.o mod_kinds.o
    306 lpm.o: modules.o cpulog.o mod_kinds.o
    307 lpm_advec.o: modules.o mod_kinds.o
    308 lpm_boundary_conds.o: modules.o cpulog.o mod_kinds.o
    309 lpm_calc_liquid_water_content.o: modules.o cpulog.o mod_kinds.o
    310 lpm_collision_kernels.o: modules.o cpulog.o user_module.o mod_kinds.o
    311 lpm_data_output_particles.o: modules.o cpulog.o mod_kinds.o
    312 lpm_droplet_collision.o: modules.o cpulog.o lpm_collision_kernels.o mod_kinds.o
    313 lpm_droplet_condensation.o: modules.o cpulog.o lpm_collision_kernels.o mod_kinds.o
    314 lpm_exchange_horiz.o: modules.o cpulog.o mod_kinds.o
    315 lpm_extend_particle_array.o: modules.o mod_kinds.o
    316 lpm_extend_tails.o: modules.o mod_kinds.o
    317 lpm_extend_tail_array.o: modules.o mod_kinds.o
    318 lpm_init.o: modules.o lpm_collision_kernels.o mod_kinds.o random_function.o
    319 lpm_init_sgs_tke.o: modules.o mod_kinds.o
    320 lpm_pack_arrays.o: modules.o mod_kinds.o
    321 lpm_read_restart_file.o: modules.o mod_kinds.o
    322 lpm_release_set.o: modules.o mod_kinds.o random_function.o
    323 lpm_set_attributes.o: modules.o cpulog.o mod_kinds.o
    324 lpm_sort_arrays.o: modules.o cpulog.o mod_kinds.o
    325 lpm_write_exchange_statistics.o: modules.o mod_kinds.o
    326 lpm_write_restart_file.o: modules.o mod_kinds.o
     308lpm.o: modules.o cpulog.o lpm_exchange_horiz.o lpm_pack_arrays.o mod_kinds.o \
     309        mod_particle_attributes.o
     310lpm_advec.o: modules.o mod_kinds.o mod_particle_attributes.o
     311lpm_boundary_conds.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
     312lpm_calc_liquid_water_content.o: modules.o cpulog.o mod_kinds.o \
     313        mod_particle_attributes.o
     314lpm_collision_kernels.o: modules.o cpulog.o user_module.o mod_kinds.o \
     315        mod_particle_attributes.o
     316lpm_data_output_particles.o: modules.o cpulog.o mod_kinds.o \
     317        mod_particle_attributes.o
     318lpm_droplet_collision.o: modules.o cpulog.o lpm_collision_kernels.o \
     319        mod_kinds.o mod_particle_attributes.o
     320lpm_droplet_condensation.o: modules.o cpulog.o lpm_collision_kernels.o \
     321        mod_kinds.o mod_particle_attributes.o
     322lpm_exchange_horiz.o: modules.o cpulog.o lpm_pack_arrays.o mod_kinds.o \
     323        mod_particle_attributes.o
     324lpm_extend_tails.o: modules.o mod_kinds.o mod_particle_attributes.o
     325lpm_extend_tail_array.o: modules.o mod_kinds.o mod_particle_attributes.o
     326lpm_init.o: modules.o lpm_collision_kernels.o mod_kinds.o \
     327        random_function.o mod_particle_attributes.o lpm_exchange_horiz.o \
     328        lpm_pack_arrays.o
     329lpm_init_sgs_tke.o: modules.o mod_kinds.o mod_particle_attributes.o
     330lpm_pack_arrays.o: modules.o mod_kinds.o mod_particle_attributes.o
     331lpm_read_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o \
     332        lpm_pack_arrays.o
     333lpm_release_set.o: modules.o mod_kinds.o random_function.o \
     334        mod_particle_attributes.o lpm_init.o
     335lpm_set_attributes.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
     336lpm_write_exchange_statistics.o: modules.o mod_kinds.o mod_particle_attributes.o
     337lpm_write_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o
    327338ls_forcing.o: modules.o cpulog.o mod_kinds.o
    328339message.o: modules.o mod_kinds.o
     
    330341modules.o: modules.f90 mod_kinds.o
    331342mod_kinds.o: mod_kinds.f90
     343mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
    332344netcdf.o: modules.o mod_kinds.o
    333345nudging.o: modules.o buoyancy.o cpulog.o mod_kinds.o
     
    343355production_e.o: modules.o mod_kinds.o wall_fluxes.o
    344356prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_u_pw.o \
    345         advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o  \
    346    advec_ws.o buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \
     357        advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o \
     358        advec_ws.o buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \
    347359        cpulog.o diffusion_e.o diffusion_s.o diffusion_u.o diffusion_v.o diffusion_w.o \
    348360        eqn_state_seawater.o impact_of_latent_heat.o mod_kinds.o microphysics.o \
    349         nudging.o plant_canopy_model.o production_e.o subsidence.o user_actions.o
     361        nudging.o plant_canopy_model.o production_e.o subsidence.o user_actions.o
    350362random_function.o: mod_kinds.o
    351363random_gauss.o: mod_kinds.o random_function.o
  • palm/trunk/SOURCE/Makefile_check

    r1321 r1359  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# mod_particle_attributes added
    2323#
    2424# Former revisions:
     
    7575      exchange_horiz_2d.f90 fft_xy.f90 init_grid.f90 init_masks.f90 \
    7676      init_cloud_physics.f90 init_pegrid.f90 local_flush.f90 local_stop.f90 \
    77       local_system.f90 message.f90 modules.f90 mod_kinds.f90 package_parin.f90 \
     77      local_system.f90 message.f90 modules.f90 mod_kinds.f90 \
     78      mod_particle_attributes.f90 package_parin.f90 \
    7879      parin.f90 poisfft.f90 random_function.f90 singleton.f90 \
    7980      subsidence.f90 temperton_fft.f90 tridia_solver.f90 \
     
    161162user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o
    162163user_last_actions.o: modules.o mod_kinds.o user_module.o
    163 user_lpm_advec.o: modules.o mod_kinds.o user_module.o
    164 user_lpm_init.o: modules.o mod_kinds.o user_module.o
    165 user_lpm_set_attributes.o: modules.o mod_kinds.o user_module.o
     164user_lpm_advec.o: modules.o mod_kinds.o mod_particle_attributes.o user_module.o
     165user_lpm_init.o: modules.o mod_kinds.o mod_particle_attributes.o user_module.o
     166user_lpm_set_attributes.o: modules.o mod_kinds.o mod_particle_attributes.o \
     167      user_module.o
    166168user_module.o: mod_kinds.o user_module.f90
    167169user_parin.o: modules.o mod_kinds.o user_module.o
  • palm/trunk/SOURCE/check_open.f90

    r1354 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Format of particle exchange statistics extended to reasonable numbers of     
     23! particles.
    2324!
    2425! Former revisions:
     
    9192
    9293    USE indices,                                                               &
    93         ONLY:  nbgp, nx, nxlg, nxrg, ny, nyng, nysg, nz, nzb
     94        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
     95               nzb, nzt
    9496
    9597    USE kinds
     
    591593!--                     unit 85 is changed (see also in routine
    592594!--                     lpm_data_output_particles)
    593              rtext = 'data format version 3.0'
     595             rtext = 'data format version 3.1'
    594596             WRITE ( 85 )  rtext
    595597             WRITE ( 85 )  number_of_particle_groups,                          &
    596598                           max_number_of_particle_groups
    597599             WRITE ( 85 )  particle_groups
     600             WRITE ( 85 )  nxl, nxr, nys, nyn, nzb, nzt, nbgp
    598601          ENDIF
    599602
     
    11361139             '#18 pt(zp)'/'#19 splptx'/'#20 splpty'/'#21 splptz')
    113711408000 FORMAT (A/                                                                &
    1138              '  step    time  # of parts   lPE sent/recv  rPE sent/recv  ',    &
    1139              'sPE sent/recv  nPE sent/recv  max # of parts'/                   &
    1140              103('-'))
     1141             '  step    time    # of parts     lPE sent/recv  rPE sent/recv  ',&
     1142             'sPE sent/recv  nPE sent/recv    max # of parts  '/               &
     1143             109('-'))
    11411144
    11421145 END SUBROUTINE check_open
  • palm/trunk/SOURCE/check_parameters.f90

    r1354 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Do not allow the execution of PALM with use_particle_tails, since particle
     23! tails are currently not supported by our new particle structure.
     24!
     25! PA0084 not necessary for new particle structure
    2326!
    2427! Former revisions:
     
    512515                        'with particle advection.'
    513516       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
     517    ENDIF
     518
     519!
     520!--
     521    IF ( use_particle_tails )  THEN
     522       message_string = 'Particle tails are currently not available due ' //   &
     523                        'to the new particle structure.'
     524       CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
    514525    ENDIF
    515526
     
    18171828
    18181829!
    1819 !-- Check the interval for sorting particles.
    1820 !-- Using particles as cloud droplets requires sorting after each timestep.
    1821     IF ( dt_sort_particles /= 0.0_wp  .AND.  cloud_droplets )  THEN
    1822        dt_sort_particles = 0.0_wp
    1823        message_string = 'dt_sort_particles is reset to 0.0 because of cloud' //&
    1824                         '_droplets = .TRUE.'
    1825        CALL message( 'check_parameters', 'PA0084', 0, 1, 0, 6, 0 )
    1826     ENDIF
    1827 
    1828 !
    18291830!-- Set the default intervals for data output, if necessary
    18301831!-- NOTE: dt_dosp has already been set in package_parin
  • palm/trunk/SOURCE/data_output_2d.f90

    r1354 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! New particle structure integrated.
    2323!
    2424! Former revisions:
     
    128128
    129129    USE particle_attributes,                                                   &
    130         ONLY:  particle_advection_start, particles, prt_count,                 &
    131                prt_start_index
     130        ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
     131               particles, prt_count
    132132   
    133133    USE pegrid
     
    163163    LOGICAL ::  two_d          !:
    164164   
    165     REAL(wp) ::  mean_r         !:
    166     REAL(wp) ::  s_r3           !:
    167     REAL(wp) ::  s_r4           !:
     165    REAL(wp) ::  mean_r        !:
     166    REAL(wp) ::  s_r2          !:
     167    REAL(wp) ::  s_r3          !:
    168168   
    169     REAL(wp), DIMENSION(:), ALLOCATABLE ::      level_z     !:
    170     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::    local_2d    !:
    171     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::    local_2d_l  !:
    172     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf    !:
     169    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  level_z             !:
     170    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  local_2d            !:
     171    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  local_2d_l          !:
     172    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf            !:
    173173    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_2d_sections   !:
    174174    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_2d_sections_l !:
     175
    175176#if defined( __parallel )
    176     REAL(wp), DIMENSION(:,:),   ALLOCATABLE ::  total_2d  !:
     177    REAL(wp), DIMENSION(:,:),   ALLOCATABLE ::  total_2d    !:
    177178#endif
    178179    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !:
     
    418419                ENDIF
    419420
    420              CASE ( 'pr_xy', 'pr_xz', 'pr_yz' )  ! mean particle radius
     421             CASE ( 'pr_xy', 'pr_xz', 'pr_yz' )  ! mean particle radius (effective radius)
    421422                IF ( av == 0 )  THEN
    422423                   IF ( simulated_time >= particle_advection_start )  THEN
     
    424425                         DO  j = nys, nyn
    425426                            DO  k = nzb, nzt+1
    426                                psi = prt_start_index(k,j,i)
     427                               number_of_particles = prt_count(k,j,i)
     428                               IF (number_of_particles <= 0)  CYCLE
     429                               particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     430                               s_r2 = 0.0_wp
    427431                               s_r3 = 0.0_wp
    428                                s_r4 = 0.0_wp
    429                                DO  n = psi, psi+prt_count(k,j,i)-1
    430                                   s_r3 = s_r3 + particles(n)%radius**3 *       &
    431                                                 particles(n)%weight_factor
    432                                   s_r4 = s_r4 + particles(n)%radius**4 *       &
    433                                                 particles(n)%weight_factor
     432                               DO  n = 1, number_of_particles
     433                                  IF ( particles(n)%particle_mask )  THEN
     434                                     s_r2 = s_r2 + particles(n)%radius**2 * &
     435                                            particles(n)%weight_factor
     436                                     s_r3 = s_r3 + particles(n)%radius**3 * &
     437                                            particles(n)%weight_factor
     438                                  ENDIF
    434439                               ENDDO
    435                                IF ( s_r3 /= 0.0_wp )  THEN
    436                                   mean_r = s_r4 / s_r3
     440                               IF ( s_r2 > 0.0_wp )  THEN
     441                                  mean_r = s_r3 / s_r2
    437442                               ELSE
    438443                                  mean_r = 0.0_wp
     
    445450                   ELSE
    446451                      tend = 0.0_wp
    447                    END IF
     452                   ENDIF
    448453                   DO  i = nxlg, nxrg
    449454                      DO  j = nysg, nyng
     
    600605                         DO  j = nys, nyn
    601606                            DO  k = nzb, nzt+1
    602                                psi = prt_start_index(k,j,i)
    603                                DO  n = psi, psi+prt_count(k,j,i)-1
    604                                   tend(k,j,i) =  tend(k,j,i) +                 &
    605                                                  particles(n)%weight_factor /  &
    606                                                  prt_count(k,j,i)
     607                               number_of_particles = prt_count(k,j,i)
     608                               IF (number_of_particles <= 0)  CYCLE
     609                               particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     610                               DO  n = 1, number_of_particles
     611                                  IF ( particles(n)%particle_mask )  THEN
     612                                     tend(k,j,i) =  tend(k,j,i) +                 &
     613                                                    particles(n)%weight_factor /  &
     614                                                    prt_count(k,j,i)
     615                                  ENDIF
    607616                               ENDDO
    608617                            ENDDO
     
    612621                   ELSE
    613622                      tend = 0.0_wp
    614                    END IF
     623                   ENDIF
    615624                   DO  i = nxlg, nxrg
    616625                      DO  j = nysg, nyng
  • palm/trunk/SOURCE/data_output_3d.f90

    r1354 r1359  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! New particle structure integrated.
    2323!
    2424! Former revisions:
     
    112112       
    113113    USE particle_attributes,                                                   &
    114         ONLY:  particles, prt_count, prt_start_index
     114        ONLY:  grid_particles, number_of_particles, particles,                 &
     115               particle_advection_start, prt_count
    115116       
    116117    USE pegrid
     
    134135
    135136    REAL(wp)     ::  mean_r    !:
     137    REAL(wp)     ::  s_r2      !:
    136138    REAL(wp)     ::  s_r3      !:
    137     REAL(wp)     ::  s_r4      !:
    138139
    139140    REAL(sp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf  !:
     
    239240          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
    240241             IF ( av == 0 )  THEN
    241                 tend = prt_count
    242                 CALL exchange_horiz( tend, nbgp )
     242                IF ( simulated_time >= particle_advection_start )  THEN
     243                   tend = prt_count
     244                   CALL exchange_horiz( tend, nbgp )
     245                ELSE
     246                   tend = 0.0_wp
     247                ENDIF
    243248                DO  i = nxlg, nxrg
    244249                   DO  j = nysg, nyng
     
    254259             ENDIF
    255260
    256           CASE ( 'pr' )  ! mean particle radius
    257              IF ( av == 0 )  THEN
    258                 DO  i = nxl, nxr
    259                    DO  j = nys, nyn
    260                       DO  k = nzb, nz_do3d
    261                          psi = prt_start_index(k,j,i)
    262                          s_r3 = 0.0_wp
    263                          s_r4 = 0.0_wp
    264                          DO  n = psi, psi+prt_count(k,j,i)-1
    265                          s_r3 = s_r3 + particles(n)%radius**3 *                &
    266                                        particles(n)%weight_factor
    267                          s_r4 = s_r4 + particles(n)%radius**4 *                &
    268                                        particles(n)%weight_factor
     261          CASE ( 'pr' )  ! mean particle radius (effective radius)
     262             IF ( av == 0 )  THEN
     263                IF ( simulated_time >= particle_advection_start )  THEN
     264                   DO  i = nxl, nxr
     265                      DO  j = nys, nyn
     266                         DO  k = nzb, nz_do3d
     267                            number_of_particles = prt_count(k,j,i)
     268                            IF (number_of_particles <= 0)  CYCLE
     269                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     270                            s_r2 = 0.0_wp
     271                            s_r3 = 0.0_wp
     272                            DO  n = 1, number_of_particles
     273                               IF ( particles(n)%particle_mask )  THEN
     274                                  s_r2 = s_r2 + particles(n)%radius**2 * &
     275                                         particles(n)%weight_factor
     276                                  s_r3 = s_r3 + particles(n)%radius**3 * &
     277                                         particles(n)%weight_factor
     278                               ENDIF
     279                            ENDDO
     280                            IF ( s_r2 > 0.0_wp )  THEN
     281                               mean_r = s_r3 / s_r2
     282                            ELSE
     283                               mean_r = 0.0_wp
     284                            ENDIF
     285                            tend(k,j,i) = mean_r
    269286                         ENDDO
    270                          IF ( s_r3 /= 0.0_wp )  THEN
    271                             mean_r = s_r4 / s_r3
    272                          ELSE
    273                             mean_r = 0.0_wp
    274                          ENDIF
    275                          tend(k,j,i) = mean_r
    276                       ENDDO
    277                    ENDDO
    278                 ENDDO
    279                 CALL exchange_horiz( tend, nbgp )
     287                      ENDDO
     288                   ENDDO
     289                   CALL exchange_horiz( tend, nbgp )
     290                ELSE
     291                   tend = 0.0_wp
     292                ENDIF
    280293                DO  i = nxlg, nxrg
    281294                   DO  j = nysg, nyng
     
    370383          CASE ( 'ql_vp' )
    371384             IF ( av == 0 )  THEN
    372                 DO  i = nxl, nxr
    373                    DO  j = nys, nyn
    374                       DO  k = nzb, nz_do3d
    375                          psi = prt_start_index(k,j,i)
    376                          DO  n = psi, psi+prt_count(k,j,i)-1
    377                             tend(k,j,i) = tend(k,j,i) +                        &
    378                                           particles(n)%weight_factor /         &
    379                                           prt_count(k,j,i)
     385                IF ( simulated_time >= particle_advection_start )  THEN
     386                   DO  i = nxl, nxr
     387                      DO  j = nys, nyn
     388                         DO  k = nzb, nz_do3d
     389                            number_of_particles = prt_count(k,j,i)
     390                            IF (number_of_particles <= 0)  CYCLE
     391                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     392                            DO  n = 1, number_of_particles
     393                               IF ( particles(n)%particle_mask )  THEN
     394                                  tend(k,j,i) =  tend(k,j,i) +                 &
     395                                                 particles(n)%weight_factor /  &
     396                                                 prt_count(k,j,i)
     397                               ENDIF
     398                            ENDDO
    380399                         ENDDO
    381400                      ENDDO
    382401                   ENDDO
    383                 ENDDO
    384                 CALL exchange_horiz( tend, nbgp )
     402                   CALL exchange_horiz( tend, nbgp )
     403                ELSE
     404                   tend = 0.0_wp
     405                ENDIF
    385406                DO  i = nxlg, nxrg
    386407                   DO  j = nysg, nyng
  • palm/trunk/SOURCE/data_output_mask.f90

    r1354 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! New particle structure integrated.
    2323!
    2424! Former revisions:
     
    9696   
    9797    USE particle_attributes,                                                   &
    98         ONLY:  particles, prt_count, prt_start_index
     98        ONLY:  grid_particles, number_of_particles, particles,                 &
     99               particle_advection_start, prt_count
    99100   
    100101    USE pegrid
     
    117118   
    118119    REAL(wp) ::  mean_r       !:
     120    REAL(wp) ::  s_r2         !:
    119121    REAL(wp) ::  s_r3         !:
    120     REAL(wp) ::  s_r4         !:
    121122   
    122123    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf    !:
     
    215216             ENDIF
    216217
    217           CASE ( 'pr' )  ! mean particle radius
    218              IF ( av == 0 )  THEN
    219                 DO  i = nxl, nxr
    220                    DO  j = nys, nyn
    221                       DO  k = nzb, nzt+1
    222                          psi = prt_start_index(k,j,i)
    223                          s_r3 = 0.0_wp
    224                          s_r4 = 0.0_wp
    225                          DO  n = psi, psi+prt_count(k,j,i)-1
    226                             s_r3 = s_r3 + particles(n)%radius**3 * &
    227                                           particles(n)%weight_factor
    228                             s_r4 = s_r4 + particles(n)%radius**4 * &
    229                                           particles(n)%weight_factor
     218          CASE ( 'pr' )  ! mean particle radius (effective radius)
     219             IF ( av == 0 )  THEN
     220                IF ( simulated_time >= particle_advection_start )  THEN
     221                   DO  i = nxl, nxr
     222                      DO  j = nys, nyn
     223                         DO  k = nzb, nz_do3d
     224                            number_of_particles = prt_count(k,j,i)
     225                            IF (number_of_particles <= 0)  CYCLE
     226                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     227                            s_r2 = 0.0_wp
     228                            s_r3 = 0.0_wp
     229                            DO  n = 1, number_of_particles
     230                               IF ( particles(n)%particle_mask )  THEN
     231                                  s_r2 = s_r2 + grid_particles(k,j,i)%particles(n)%radius**2 * &
     232                                         grid_particles(k,j,i)%particles(n)%weight_factor
     233                                  s_r3 = s_r3 + grid_particles(k,j,i)%particles(n)%radius**3 * &
     234                                         grid_particles(k,j,i)%particles(n)%weight_factor
     235                               ENDIF
     236                            ENDDO
     237                            IF ( s_r2 > 0.0_wp )  THEN
     238                               mean_r = s_r3 / s_r2
     239                            ELSE
     240                               mean_r = 0.0_wp
     241                            ENDIF
     242                            tend(k,j,i) = mean_r
    230243                         ENDDO
    231                          IF ( s_r3 /= 0.0_wp )  THEN
    232                             mean_r = s_r4 / s_r3
    233                          ELSE
    234                             mean_r = 0.0_wp
    235                          ENDIF
    236                          tend(k,j,i) = mean_r
    237                       ENDDO
    238                    ENDDO
    239                 ENDDO
    240                 CALL exchange_horiz( tend, nbgp )
     244                      ENDDO
     245                   ENDDO
     246                   CALL exchange_horiz( tend, nbgp )
     247                ELSE
     248                   tend = 0.0_wp
     249                ENDIF
    241250                DO  i = 1, mask_size_l(mid,1)
    242251                   DO  j = 1, mask_size_l(mid,2)
     
    304313          CASE ( 'ql_vp' )
    305314             IF ( av == 0 )  THEN
    306                 DO  i = nxl, nxr
    307                    DO  j = nys, nyn
    308                       DO  k = nzb, nz_do3d
    309                          psi = prt_start_index(k,j,i)
    310                          DO  n = psi, psi+prt_count(k,j,i)-1
    311                             tend(k,j,i) = tend(k,j,i) + &
     315                IF ( simulated_time >= particle_advection_start )  THEN
     316                   DO  i = nxl, nxr
     317                      DO  j = nys, nyn
     318                         DO  k = nzb, nz_do3d
     319                            number_of_particles = prt_count(k,j,i)
     320                            IF (number_of_particles <= 0)  CYCLE
     321                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     322                            DO  n = 1, number_of_particles
     323                               IF ( particles(n)%particle_mask )  THEN
     324                                  tend(k,j,i) = tend(k,j,i) + &
    312325                                          particles(n)%weight_factor / &
    313326                                          prt_count(k,j,i)
     327                               ENDIF
     328                            ENDDO
    314329                         ENDDO
    315330                      ENDDO
    316331                   ENDDO
    317                 ENDDO
    318                 CALL exchange_horiz( tend, nbgp )
     332                   CALL exchange_horiz( tend, nbgp )
     333                ELSE
     334                   tend = 0.0_wp
     335                ENDIF
    319336                DO  i = 1, mask_size_l(mid,1)
    320337                   DO  j = 1, mask_size_l(mid,2)
  • palm/trunk/SOURCE/data_output_ptseries.f90

    r1354 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! New particle structure integrated.
    2323!
    2424! Former revisions:
     
    7070
    7171    USE indices,                                                               &
    72         ONLY:
     72        ONLY: nxl, nxr, nys, nyn, nzb, nzt
    7373
    7474    USE kinds
     
    7777
    7878    USE particle_attributes,                                                   &
    79         ONLY:  number_of_particles, number_of_particle_groups, particles
     79        ONLY:  grid_particles, number_of_particles, number_of_particle_groups, &
     80               particles, prt_count
    8081
    8182    USE pegrid
     
    8788    INTEGER(iwp) ::  inum !:
    8889    INTEGER(iwp) ::  j    !:
     90    INTEGER(iwp) ::  jg   !:
     91    INTEGER(iwp) ::  k    !:
    8992    INTEGER(iwp) ::  n    !:
    9093
     
    120123!-- Calculate or collect the particle time series quantities for all particles
    121124!-- and seperately for each particle group (if there is more than one group)
    122     DO  n = 1, number_of_particles
    123 
    124        pts_value_l(0,1)  = number_of_particles            ! total # of particles
    125        pts_value_l(0,2)  = pts_value_l(0,2) + &
    126                            ( particles(n)%x - particles(n)%origin_x )  ! mean x
    127        pts_value_l(0,3)  = pts_value_l(0,3) + &
    128                            ( particles(n)%y - particles(n)%origin_y )  ! mean y
    129        pts_value_l(0,4)  = pts_value_l(0,4) + &
    130                            ( particles(n)%z - particles(n)%origin_z )  ! mean z
    131        pts_value_l(0,5)  = pts_value_l(0,5) + particles(n)%z ! mean z (absolute)
    132        pts_value_l(0,6)  = pts_value_l(0,6) + particles(n)%speed_x     ! mean u
    133        pts_value_l(0,7)  = pts_value_l(0,7) + particles(n)%speed_y     ! mean v
    134        pts_value_l(0,8)  = pts_value_l(0,8) + particles(n)%speed_z     ! mean w
    135        IF ( .NOT. curvature_solution_effects )  THEN
    136           pts_value_l(0,9)  = pts_value_l(0,9) + particles(n)%rvar1  ! mean sgsu
    137           pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
    138           pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
    139        ENDIF
    140        IF ( particles(n)%speed_z > 0.0_wp )  THEN
    141           pts_value_l(0,12) = pts_value_l(0,12) + 1.0_wp  ! # of upward moving prts
    142           pts_value_l(0,13) = pts_value_l(0,13) + &
    143                                             particles(n)%speed_z ! mean w upw.
    144        ELSE
    145           pts_value_l(0,14) = pts_value_l(0,14) + &
    146                                             particles(n)%speed_z ! mean w down
    147        ENDIF
    148        pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad
    149        pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad
    150        pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad
    151        pts_value_l(0,18) = number_of_particles
    152        pts_value_l(0,19) = number_of_particles
    153 
    154 !
    155 !--    Repeat the same for the respective particle group
    156        IF ( number_of_particle_groups > 1 )  THEN
    157           j = particles(n)%group
    158 
    159           pts_value_l(j,1)  = pts_value_l(j,1)  + 1
    160           pts_value_l(j,2)  = pts_value_l(j,2)  + &
    161                               ( particles(n)%x - particles(n)%origin_x )
    162           pts_value_l(j,3)  = pts_value_l(j,3)  + &
    163                               ( particles(n)%y - particles(n)%origin_y )
    164           pts_value_l(j,4)  = pts_value_l(j,4)  + &
    165                               ( particles(n)%z - particles(n)%origin_z )
    166           pts_value_l(j,5)  = pts_value_l(j,5)  + particles(n)%z
    167           pts_value_l(j,6)  = pts_value_l(j,6)  + particles(n)%speed_x
    168           pts_value_l(j,7)  = pts_value_l(j,7)  + particles(n)%speed_y
    169           pts_value_l(j,8)  = pts_value_l(j,8)  + particles(n)%speed_z
    170           IF ( .NOT. curvature_solution_effects )  THEN
    171              pts_value_l(j,9)  = pts_value_l(j,9)  + particles(n)%rvar1
    172              pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2
    173              pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3
    174           ENDIF
    175           IF ( particles(n)%speed_z > 0.0_wp )  THEN
    176              pts_value_l(j,12) = pts_value_l(j,12) + 1.0_wp
    177              pts_value_l(j,13) = pts_value_l(j,13) + particles(n)%speed_z
    178           ELSE
    179              pts_value_l(j,14) = pts_value_l(j,14) + particles(n)%speed_z
    180           ENDIF
    181           pts_value_l(j,15) = pts_value_l(j,15) + particles(n)%radius
    182           pts_value_l(j,16) = MIN( pts_value(j,16), particles(n)%radius )
    183           pts_value_l(j,17) = MAX( pts_value(j,17), particles(n)%radius )
    184           pts_value_l(j,18) = pts_value_l(j,18) + 1.0_wp
    185           pts_value_l(j,19) = pts_value_l(j,19) + 1.0_wp
    186 
    187        ENDIF
    188 
     125    DO  i = nxl, nxr
     126       DO  j = nys, nyn
     127          DO  k = nzb, nzt
     128             number_of_particles = prt_count(k,j,i)
     129             IF (number_of_particles <= 0)  CYCLE
     130             particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     131             DO  n = 1, number_of_particles
     132
     133                IF ( particles(n)%particle_mask )  THEN  ! Restrict analysis to active particles
     134
     135                   pts_value_l(0,1)  = pts_value_l(0,1) + 1.0_wp  ! total # of particles
     136                   pts_value_l(0,2)  = pts_value_l(0,2) +                      &
     137                          ( particles(n)%x - particles(n)%origin_x )  ! mean x
     138                   pts_value_l(0,3)  = pts_value_l(0,3) +                      &
     139                          ( particles(n)%y - particles(n)%origin_y )  ! mean y
     140                   pts_value_l(0,4)  = pts_value_l(0,4) +                      &
     141                          ( particles(n)%z - particles(n)%origin_z )  ! mean z
     142                   pts_value_l(0,5)  = pts_value_l(0,5) + particles(n)%z        ! mean z (absolute)
     143                   pts_value_l(0,6)  = pts_value_l(0,6) + particles(n)%speed_x  ! mean u
     144                   pts_value_l(0,7)  = pts_value_l(0,7) + particles(n)%speed_y  ! mean v
     145                   pts_value_l(0,8)  = pts_value_l(0,8) + particles(n)%speed_z  ! mean w
     146                   IF ( .NOT. curvature_solution_effects )  THEN
     147                      pts_value_l(0,9)  = pts_value_l(0,9)  + particles(n)%rvar1 ! mean sgsu
     148                      pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
     149                      pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
     150                   ENDIF
     151                   IF ( particles(n)%speed_z > 0.0_wp )  THEN
     152                      pts_value_l(0,12) = pts_value_l(0,12) + 1.0_wp  ! # of upward moving prts
     153                      pts_value_l(0,13) = pts_value_l(0,13) +                  &
     154                                              particles(n)%speed_z ! mean w upw.
     155                   ELSE
     156                      pts_value_l(0,14) = pts_value_l(0,14) +                  &
     157                                              particles(n)%speed_z ! mean w down
     158                   ENDIF
     159                   pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad
     160                   pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad
     161                   pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad
     162                   pts_value_l(0,18) = pts_value_l(0,18) + 1.0_wp
     163                   pts_value_l(0,19) = pts_value_l(0,18) + 1.0_wp
     164!
     165!--                Repeat the same for the respective particle group
     166                   IF ( number_of_particle_groups > 1 )  THEN
     167                      jg = particles(n)%group
     168
     169                      pts_value_l(jg,1)  = pts_value_l(jg,1) + 1.0_wp
     170                      pts_value_l(jg,2)  = pts_value_l(jg,2) +                   &
     171                           ( particles(n)%x - particles(n)%origin_x )
     172                      pts_value_l(jg,3)  = pts_value_l(jg,3) +                   &
     173                           ( particles(n)%y - particles(n)%origin_y )
     174                      pts_value_l(jg,4)  = pts_value_l(jg,4) +                   &
     175                           ( particles(n)%z - particles(n)%origin_z )
     176                      pts_value_l(jg,5)  = pts_value_l(jg,5) + particles(n)%z
     177                      pts_value_l(jg,6)  = pts_value_l(jg,6) + particles(n)%speed_x
     178                      pts_value_l(jg,7)  = pts_value_l(jg,7) + particles(n)%speed_y
     179                      pts_value_l(jg,8)  = pts_value_l(jg,8) + particles(n)%speed_z
     180                      IF ( .NOT. curvature_solution_effects )  THEN
     181                         pts_value_l(jg,9)  = pts_value_l(jg,9)  + particles(n)%rvar1
     182                         pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
     183                         pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
     184                      ENDIF
     185                      IF ( particles(n)%speed_z > 0.0_wp )  THEN
     186                         pts_value_l(jg,12) = pts_value_l(jg,12) + 1.0_wp
     187                         pts_value_l(jg,13) = pts_value_l(jg,13) + particles(n)%speed_z
     188                      ELSE
     189                         pts_value_l(jg,14) = pts_value_l(jg,14) + particles(n)%speed_z
     190                      ENDIF
     191                      pts_value_l(jg,15) = pts_value_l(jg,15) + particles(n)%radius
     192                      pts_value_l(jg,16) = MIN( pts_value(jg,16), particles(n)%radius )
     193                      pts_value_l(jg,17) = MAX( pts_value(jg,17), particles(n)%radius )
     194                      pts_value_l(jg,18) = pts_value_l(jg,18) + 1.0_wp
     195                      pts_value_l(jg,19) = pts_value_l(jg,19) + 1.0_wp
     196                   ENDIF
     197
     198                ENDIF
     199
     200             ENDDO
     201
     202          ENDDO
     203       ENDDO
    189204    ENDDO
     205
    190206
    191207#if defined( __parallel )
     
    243259!-- Calculate higher order moments of particle time series quantities,
    244260!-- seperately for each particle group (if there is more than one group)
    245     DO  n = 1, number_of_particles
    246 
    247        pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - &
    248                            particles(n)%origin_x - pts_value(0,2) )**2 ! x*2
    249        pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - &
    250                            particles(n)%origin_y - pts_value(0,3) )**2 ! y*2
    251        pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - &
    252                            particles(n)%origin_z - pts_value(0,4) )**2 ! z*2
    253        pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - &
    254                                                 pts_value(0,6) )**2   ! u*2
    255        pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - &
    256                                                  pts_value(0,7) )**2   ! v*2
    257        pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
    258                                                  pts_value(0,8) )**2   ! w*2
    259        IF ( .NOT. curvature_solution_effects )  THEN
    260           pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
    261                                                     pts_value(0,9) )**2   ! u"2
    262           pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
    263                                                     pts_value(0,10) )**2  ! v"2
    264           pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
    265                                                     pts_value(0,11) )**2  ! w"2
    266        ENDIF
    267 !
    268 !--    Repeat the same for the respective particle group
    269        IF ( number_of_particle_groups > 1 )  THEN
    270           j = particles(n)%group
    271 
    272           pts_value_l(j,20) = pts_value_l(j,20) + ( particles(n)%x - &
    273                               particles(n)%origin_x - pts_value(j,2) )**2
    274           pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%y - &
    275                               particles(n)%origin_y - pts_value(j,3) )**2
    276           pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%z - &
    277                               particles(n)%origin_z - pts_value(j,4) )**2
    278           pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_x - &
    279                                                     pts_value(j,6) )**2
    280           pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%speed_y - &
    281                                                     pts_value(j,7) )**2
    282           pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%speed_z - &
    283                                                     pts_value(j,8) )**2
    284           IF ( .NOT. curvature_solution_effects )  THEN
    285              pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar1 - &
    286                                                        pts_value(j,9) )**2
    287              pts_value_l(j,27) = pts_value_l(j,27) + ( particles(n)%rvar2 - &
    288                                                        pts_value(j,10) )**2
    289              pts_value_l(j,28) = pts_value_l(j,28) + ( particles(n)%rvar3 - &
    290                                                        pts_value(j,11) )**2
    291           ENDIF
    292        ENDIF
    293 
     261    DO  i = nxl, nxr
     262       DO  j = nys, nyn
     263          DO  k = nzb, nzt
     264             number_of_particles = prt_count(k,j,i)
     265             IF (number_of_particles <= 0)  CYCLE
     266             particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     267             DO  n = 1, number_of_particles
     268
     269                pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - &
     270                                    particles(n)%origin_x - pts_value(0,2) )**2 ! x*2
     271                pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - &
     272                                    particles(n)%origin_y - pts_value(0,3) )**2 ! y*2
     273                pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - &
     274                                    particles(n)%origin_z - pts_value(0,4) )**2 ! z*2
     275                pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - &
     276                                                         pts_value(0,6) )**2   ! u*2
     277                pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - &
     278                                                          pts_value(0,7) )**2   ! v*2
     279                pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
     280                                                          pts_value(0,8) )**2   ! w*2
     281                IF ( .NOT. curvature_solution_effects )  THEN
     282                   pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
     283                                                             pts_value(0,9) )**2   ! u"2
     284                   pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
     285                                                             pts_value(0,10) )**2  ! v"2
     286                   pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
     287                                                             pts_value(0,11) )**2  ! w"2
     288                ENDIF
     289!
     290!--             Repeat the same for the respective particle group
     291                IF ( number_of_particle_groups > 1 )  THEN
     292                   jg = particles(n)%group
     293
     294                   pts_value_l(jg,20) = pts_value_l(jg,20) + ( particles(n)%x - &
     295                                       particles(n)%origin_x - pts_value(jg,2) )**2
     296                   pts_value_l(jg,21) = pts_value_l(jg,21) + ( particles(n)%y - &
     297                                       particles(n)%origin_y - pts_value(jg,3) )**2
     298                   pts_value_l(jg,22) = pts_value_l(jg,22) + ( particles(n)%z - &
     299                                       particles(n)%origin_z - pts_value(jg,4) )**2
     300                   pts_value_l(jg,23) = pts_value_l(jg,23) + ( particles(n)%speed_x - &
     301                                                             pts_value(jg,6) )**2
     302                   pts_value_l(jg,24) = pts_value_l(jg,24) + ( particles(n)%speed_y - &
     303                                                             pts_value(jg,7) )**2
     304                   pts_value_l(jg,25) = pts_value_l(jg,25) + ( particles(n)%speed_z - &
     305                                                             pts_value(jg,8) )**2
     306                   IF ( .NOT. curvature_solution_effects )  THEN
     307                      pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
     308                                                                pts_value(jg,9) )**2
     309                      pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
     310                                                                pts_value(jg,10) )**2
     311                      pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
     312                                                                pts_value(jg,11) )**2
     313                   ENDIF
     314                ENDIF
     315
     316             ENDDO
     317          ENDDO
     318       ENDDO
    294319    ENDDO
    295320
  • palm/trunk/SOURCE/header.f90

    r1354 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! dt_sort_particles removed
    2323!
    2424! Former revisions:
     
    174174        ONLY:  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, collision_kernel,     &
    175175               density_ratio, dissipation_classes, dt_min_part, dt_prel,       &
    176                dt_sort_particles, dt_write_particle_data, end_time_prel,       &
     176               dt_write_particle_data, end_time_prel,                          &
    177177               maximum_number_of_tailpoints, maximum_tailpoint_age,            &
    178178               minimum_tailpoint_distance, number_of_particle_groups,          &
     
    15761576       WRITE ( io, 480 )  particle_advection_start, dt_prel, bc_par_lr, &
    15771577                          bc_par_ns, bc_par_b, bc_par_t, particle_maximum_age, &
    1578                           end_time_prel, dt_sort_particles
     1578                          end_time_prel
    15791579       IF ( use_sgs_for_particles )  WRITE ( io, 488 )  dt_min_part
    15801580       IF ( random_start_position )  WRITE ( io, 481 )
     
    20472047            '                            bottom:     ', A, ' top:         ', A/&
    20482048            '       Maximum particle age:                 ',F9.1,' s'/ &
    2049             '       Advection stopped at t = ',F9.1,' s'/ &
    2050             '       Particles are sorted every ',F9.1,' s'/)
     2049            '       Advection stopped at t = ',F9.1,' s'/)
    20512050481 FORMAT ('       Particles have random start positions'/)
    20522051482 FORMAT ('          Particles are advected only horizontally'/)
  • palm/trunk/SOURCE/init_3d_model.f90

    r1354 r1359  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! module lpm_init_mod added to use statements, because lpm_init has become a
     23! module
    2324!
    2425! Former revisions:
     
    168169   
    169170    USE indices
     171
     172    USE lpm_init_mod,                                                              &
     173        ONLY:  lpm_init
    170174   
    171175    USE kinds
  • palm/trunk/SOURCE/lpm.f90

    r1321 r1359  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2324!
    2425! Former revisions:
     
    8283    USE control_parameters,                                                    &
    8384        ONLY:  cloud_droplets, dt_3d, dt_3d_reached, dt_3d_reached_l,          &
    84                molecular_viscosity, simulated_time
     85               molecular_viscosity, simulated_time, topography
    8586
    8687    USE cpulog,                                                                &
    8788        ONLY:  cpu_log, log_point, log_point_s
    8889
     90    USE indices,                                                               &
     91        ONLY: nxl, nxr, nys, nyn, nzb, nzt
     92
    8993    USE kinds
    9094
     95    USE lpm_exchange_horiz_mod,                                                &
     96        ONLY:  lpm_exchange_horiz, lpm_move_particle
     97
     98    USE lpm_pack_arrays_mod,                                                   &
     99        ONLY:  lpm_pack_all_arrays
     100
    91101    USE particle_attributes,                                                   &
    92         ONLY:  collision_kernel, deleted_particles, dt_sort_particles,         &
    93                deleted_tails, dt_write_particle_data, dt_prel, end_time_prel,  &
    94                number_of_particles, number_of_particle_groups,particles,      &
    95                particle_groups, particle_mask, trlp_count_sum, tail_mask,      &
    96                time_prel, time_sort_particles, time_write_particle_data,       &
    97                trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
    98                trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
    99                trsp_count_recv_sum, use_particle_tails, use_sgs_for_particles, &
    100                write_particle_statistics
     102        ONLY:  collision_kernel, deleted_particles, deleted_tails,             &
     103               dt_write_particle_data, dt_prel, end_time_prel,                 &
     104               grid_particles, number_of_particles, number_of_particle_groups, &
     105               particles, particle_groups, prt_count, trlp_count_sum,          &
     106               tail_mask, time_prel, time_sort_particles,                      &
     107               time_write_particle_data, trlp_count_recv_sum, trnp_count_sum,  &
     108               trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &
     109               trsp_count_sum, trsp_count_recv_sum, use_particle_tails,        &
     110               use_sgs_for_particles, write_particle_statistics
    101111
    102112    USE pegrid
     
    104114    IMPLICIT NONE
    105115
    106     INTEGER(iwp) ::  m                       !:
    107 
     116    INTEGER(iwp)       ::  i                  !:
     117    INTEGER(iwp)       ::  ie                 !:
     118    INTEGER(iwp)       ::  is                 !:
     119    INTEGER(iwp)       ::  j                  !:
     120    INTEGER(iwp)       ::  je                 !:
     121    INTEGER(iwp)       ::  js                 !:
     122    INTEGER(iwp)       ::  k                  !:
     123    INTEGER(iwp)       ::  ke                 !:
     124    INTEGER(iwp)       ::  ks                 !:
     125    INTEGER(iwp)       ::  m                  !:
     126    INTEGER(iwp), SAVE ::  steps = 0          !:
     127
     128    LOGICAL            ::  first_loop_stride  !:
    108129
    109130    CALL cpu_log( log_point(25), 'lpm', 'start' )
     
    125146    ENDIF
    126147
    127 
    128148!
    129149!-- Initialize arrays for marking those particles/tails to be deleted after the
    130150!-- (sub-) timestep
    131     particle_mask     = .TRUE.
    132151    deleted_particles = 0
    133152
     
    157176!-- of the particle groups
    158177    DO  m = 1, number_of_particle_groups
    159        IF ( particle_groups(m)%density_ratio /= 0.0 )  THEN
     178       IF ( particle_groups(m)%density_ratio /= 0.0_wp )  THEN
    160179          particle_groups(m)%exp_arg  =                                        &
    161                     4.5 * particle_groups(m)%density_ratio *                   &
     180                    4.5_wp * particle_groups(m)%density_ratio *                &
    162181                    molecular_viscosity / ( particle_groups(m)%radius )**2
    163           particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
    164                                              dt_3d )
     182
     183          particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg *     &
     184                    dt_3d )
    165185       ENDIF
    166186    ENDDO
    167187
    168 
    169 !
    170 !-- Particle (droplet) growth by condensation/evaporation and collision
    171     IF ( cloud_droplets )  THEN
    172 
    173 !
    174 !--    Reset summation arrays
    175        ql_c = 0.0;  ql_v = 0.0;  ql_vp = 0.0
    176 
    177 !
    178 !--    Droplet growth by condensation / evaporation
    179        CALL lpm_droplet_condensation
    180 
    181 !
    182 !--    Particle growth by collision
    183        IF ( collision_kernel /= 'none' )  CALL lpm_droplet_collision
    184 
    185     ENDIF
    186 
    187 
    188 !
    189 !-- If particle advection includes SGS velocity components, calculate the
    190 !-- required SGS quantities (i.e. gradients of the TKE, as well as horizontally
    191 !-- averaged profiles of the SGS TKE and the resolved-scale velocity variances)
    192     IF ( use_sgs_for_particles )  CALL lpm_init_sgs_tke
    193 
    194 
    195 !
    196 !-- Initialize the variable storing the total time that a particle has advanced
    197 !-- within the timestep procedure
    198     particles(1:number_of_particles)%dt_sum = 0.0
    199 
    200 
     188!
     189!-- If necessary, release new set of particles
     190    IF ( time_prel >= dt_prel  .AND.  end_time_prel > simulated_time )  THEN
     191
     192       CALL lpm_release_set
     193!
     194!--    The MOD function allows for changes in the output interval with
     195!--    restart runs.
     196       time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
     197
     198    ENDIF
     199!
     200!-- Reset summation arrays
     201    IF ( cloud_droplets)  THEN
     202       ql_c  = 0.0_wp
     203       ql_v  = 0.0_wp
     204       ql_vp = 0.0_wp
     205    ENDIF
     206
     207    first_loop_stride = .TRUE.
     208    grid_particles(:,:,:)%time_loop_done = .TRUE.
    201209!
    202210!-- Timestep loop for particle advection.
     
    204212!-- (within the total domain!) has reached the LES timestep (dt_3d).
    205213!-- In case of including the SGS velocities, the particle timestep may be
    206 !-- smaller than the LES timestep (because of the Lagrangian timescale restric-
    207 !-- tion) and particles may require to undergo several particle timesteps,
    208 !-- before the LES timestep is reached. Because the number of these particle
    209 !-- timesteps to be carried out is unknown at first, these steps are carried
    210 !-- out in the following infinite loop with exit condition.
     214!-- smaller than the LES timestep (because of the Lagrangian timescale
     215!-- restriction) and particles may require to undergo several particle
     216!-- timesteps, before the LES timestep is reached. Because the number of these
     217!-- particle timesteps to be carried out is unknown at first, these steps are
     218!-- carried out in the following infinite loop with exit condition.
    211219    DO
    212 
    213220       CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
    214 
    215 !
    216 !--    Initialize the switch used for the loop exit condition checked at the
    217 !--    end of this loop.
    218 !--    If at least one particle has failed to reach the LES timestep, this
    219 !--    switch will be set false.
    220        dt_3d_reached_l = .TRUE.
    221 
    222 !
    223 !--    Particle advection
    224        CALL lpm_advec
    225 
    226 !
    227 !--    Particle reflection from walls
    228        CALL lpm_boundary_conds( 'walls' )
    229 
    230 !
    231 !--    User-defined actions after the calculation of the new particle position
    232        CALL user_lpm_advec
     221       CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
     222!
     223!--    If particle advection includes SGS velocity components, calculate the
     224!--    required SGS quantities (i.e. gradients of the TKE, as well as
     225!--    horizontally averaged profiles of the SGS TKE and the resolved-scale
     226!--    velocity variances)
     227
     228       IF ( use_sgs_for_particles )  CALL lpm_init_sgs_tke
     229
     230       DO  i = nxl, nxr
     231          DO  j = nys, nyn
     232             DO  k = nzb+1, nzt
     233
     234                number_of_particles = prt_count(k,j,i)
     235!
     236!--             If grid cell gets empty, flag must be true
     237                IF ( number_of_particles <= 0 )  THEN
     238                   grid_particles(k,j,i)%time_loop_done = .TRUE.
     239                   CYCLE
     240                ENDIF
     241
     242                IF ( .NOT. first_loop_stride  .AND.  &
     243                     grid_particles(k,j,i)%time_loop_done ) CYCLE
     244
     245                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     246
     247                particles(1:number_of_particles)%particle_mask = .TRUE.
     248!
     249!--             Initialize the variable storing the total time that a particle
     250!--             has advanced within the timestep procedure
     251                IF ( first_loop_stride )  THEN
     252                   particles(1:number_of_particles)%dt_sum = 0.0_wp
     253                ENDIF
     254!
     255!--             Particle (droplet) growth by condensation/evaporation and
     256!--             collision
     257                IF ( cloud_droplets  .AND.  first_loop_stride)  THEN
     258!
     259!--                Droplet growth by condensation / evaporation
     260                   CALL lpm_droplet_condensation(i,j,k)
     261!
     262!--                Particle growth by collision
     263                   IF ( collision_kernel /= 'none' )  THEN
     264                      CALL lpm_droplet_collision(i,j,k)
     265                   ENDIF
     266
     267                ENDIF
     268!
     269!--             Initialize the switch used for the loop exit condition checked
     270!--             at the end of this loop. If at least one particle has failed to
     271!--             reach the LES timestep, this switch will be set false in
     272!--             lpm_advec.
     273                dt_3d_reached_l = .TRUE.
     274
     275!
     276!--             Particle advection
     277                CALL lpm_advec(i,j,k)
     278!
     279!--             Particle reflection from walls
     280                IF ( topography /= 'flat' )  THEN
     281                   CALL lpm_boundary_conds( 'walls' )
     282                ENDIF
     283!
     284!--             User-defined actions after the calculation of the new particle
     285!--             position
     286                CALL user_lpm_advec
     287!
     288!--             Apply boundary conditions to those particles that have crossed
     289!--             the top or bottom boundary and delete those particles, which are
     290!--             older than allowed
     291                CALL lpm_boundary_conds( 'bottom/top' )
     292!
     293!---            If not all particles of the actual grid cell have reached the
     294!--             LES timestep, this cell has to to another loop iteration. Due to
     295!--             the fact that particles can move into neighboring grid cell,
     296!--             these neighbor cells also have to perform another loop iteration
     297                IF ( .NOT. dt_3d_reached_l )  THEN
     298                   ks = MAX(nzb+1,k)
     299                   ke = MIN(nzt,k)
     300                   js = MAX(nys,j)
     301                   je = MIN(nyn,j)
     302                   is = MAX(nxl,i)
     303                   ie = MIN(nxr,i)
     304                   grid_particles(ks:ke,js:je,is:ie)%time_loop_done = .FALSE.
     305                ENDIF
     306
     307             ENDDO
     308          ENDDO
     309       ENDDO
     310
     311       steps = steps + 1
     312       dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
    233313
    234314!
     
    250330
    251331!
    252 !--    If necessary, release new set of particles
    253        IF ( time_prel >= dt_prel  .AND.  end_time_prel > simulated_time  .AND. &
    254             dt_3d_reached )  THEN
    255 
    256           CALL lpm_release_set
    257 
    258 !
    259 !--       The MOD function allows for changes in the output interval with
    260 !--       restart runs.
    261           time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
    262 
    263        ENDIF
     332!--    Move Particles local to PE to a different grid cell
     333       CALL lpm_move_particle
    264334
    265335!
     
    268338
    269339!
    270 !--    Apply boundary conditions to those particles that have crossed the top or
    271 !--    bottom boundary and delete those particles, which are older than allowed
    272        CALL lpm_boundary_conds( 'bottom/top' )
    273 
    274 !
    275340!--    Pack particles (eliminate those marked for deletion),
    276341!--    determine new number of particles
    277        IF ( number_of_particles > 0  .AND.  deleted_particles > 0 )  THEN
    278           CALL lpm_pack_arrays
    279        ENDIF
    280 
    281 !
    282 !--    Initialize variables for the next (sub-) timestep, i.e. for marking those
    283 !--    particles to be deleted after the timestep
    284        particle_mask     = .TRUE.
     342       CALL lpm_pack_all_arrays
     343
     344!
     345!--    Initialize variables for the next (sub-) timestep, i.e., for marking
     346!--    those particles to be deleted after the timestep
    285347       deleted_particles = 0
    286348
     
    293355       IF ( dt_3d_reached )  EXIT
    294356
     357       first_loop_stride = .FALSE.
    295358    ENDDO   ! timestep loop
    296359
    297 
    298 !
    299 !-- Sort particles in the sequence the gridboxes are stored in the memory
    300     time_sort_particles = time_sort_particles + dt_3d
    301     IF ( time_sort_particles >= dt_sort_particles )  THEN
    302        CALL lpm_sort_arrays
    303        time_sort_particles = MOD( time_sort_particles, &
    304                                   MAX( dt_sort_particles, dt_3d ) )
    305     ENDIF
    306 
    307 
    308360!
    309361!-- Calculate the new liquid water content for each grid box
    310     IF ( cloud_droplets )  CALL lpm_calc_liquid_water_content
     362    IF ( cloud_droplets )  THEN
     363       CALL lpm_calc_liquid_water_content
     364    ENDIF
     365
    311366
    312367
     
    325380
    326381!
     382!-- particle tails currently not available
     383!
    327384!-- If required, add the current particle positions to the particle tails
    328     IF ( use_particle_tails )  CALL lpm_extend_tails
     385!   IF ( use_particle_tails )  CALL lpm_extend_tails
    329386
    330387
     
    336393    CALL cpu_log( log_point(25), 'lpm', 'stop' )
    337394
    338 
    339395 END SUBROUTINE lpm
  • palm/trunk/SOURCE/lpm_advec.f90

    r1323 r1359  
    1  SUBROUTINE lpm_advec
     1 SUBROUTINE lpm_advec (ip,jp,kp)
    22
    33!--------------------------------------------------------------------------------!
     
    2020! Current revisions:
    2121! ------------------
    22 !
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2324!
    2425! Former revisions:
     
    5657
    5758    USE arrays_3d,                                                             &
    58         ONLY:  de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu, zw
     59        ONLY:  de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu,  &
     60               zw
     61
     62    USE cpulog
     63
     64    USE interfaces
     65
     66    USE pegrid
    5967
    6068    USE control_parameters,                                                    &
    6169        ONLY:  atmos_ocean_sign, cloud_droplets, dt_3d, dt_3d_reached_l, dz,   &
    6270               g, kappa, molecular_viscosity, prandtl_layer, topography,       &
    63                u_gtrans, v_gtrans
     71               u_gtrans, v_gtrans, simulated_time
    6472
    6573    USE grid_variables,                                                        &
     
    7280   
    7381    USE particle_attributes,                                                   &
    74         ONLY:  c_0, density_ratio, dt_min_part, iran_part, log_z_z0,           &
    75                number_of_particles, number_of_sublayers, particles,            &
    76                particle_groups, offset_ocean_nzt, offset_ocean_nzt_m1,         &
    77                sgs_wfu_part, sgs_wfv_part, sgs_wfw_part, use_sgs_for_particles,&
    78                vertical_particle_advection, z0_av_global
     82        ONLY:  block_offset, c_0, density_ratio, dt_min_part, grid_particles,  &
     83               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
     84               particles, particle_groups, offset_ocean_nzt,                   &
     85               offset_ocean_nzt_m1, sgs_wfu_part, sgs_wfv_part, sgs_wfw_part,  &
     86               use_sgs_for_particles, vertical_particle_advection, z0_av_global
    7987       
    8088    USE statistics,                                                            &
    8189        ONLY:  hom
    82        
    8390
    8491    IMPLICIT NONE
     
    8794    INTEGER(iwp) ::  gp_outside_of_building(1:8) !:
    8895    INTEGER(iwp) ::  i                           !:
     96    INTEGER(iwp) ::  ip                          !:
    8997    INTEGER(iwp) ::  j                           !:
     98    INTEGER(iwp) ::  jp                          !:
    9099    INTEGER(iwp) ::  k                           !:
     100    INTEGER(iwp) ::  kp                          !:
    91101    INTEGER(iwp) ::  kw                          !:
    92102    INTEGER(iwp) ::  n                           !:
     103    INTEGER(iwp) ::  nb                          !:
    93104    INTEGER(iwp) ::  num_gp                      !:
     105
     106    INTEGER(iwp), DIMENSION(0:7) ::  start_index !:
     107    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !:
    94108
    95109    REAL(wp) ::  aa                 !:
     
    99113    REAL(wp) ::  d_z_p_z0           !:
    100114    REAL(wp) ::  dd                 !:
    101     REAL(wp) ::  de_dx_int          !:
    102115    REAL(wp) ::  de_dx_int_l        !:
    103116    REAL(wp) ::  de_dx_int_u        !:
    104     REAL(wp) ::  de_dy_int          !:
    105117    REAL(wp) ::  de_dy_int_l        !:
    106118    REAL(wp) ::  de_dy_int_u        !:
    107119    REAL(wp) ::  de_dt              !:
    108120    REAL(wp) ::  de_dt_min          !:
    109     REAL(wp) ::  de_dz_int          !:
    110121    REAL(wp) ::  de_dz_int_l        !:
    111122    REAL(wp) ::  de_dz_int_u        !:
    112     REAL(wp) ::  dens_ratio         !:
    113     REAL(wp) ::  diss_int           !:
    114123    REAL(wp) ::  diss_int_l         !:
    115124    REAL(wp) ::  diss_int_u         !:
    116125    REAL(wp) ::  dt_gap             !:
    117     REAL(wp) ::  dt_particle        !:
    118126    REAL(wp) ::  dt_particle_m      !:
    119     REAL(wp) ::  e_int              !:
    120127    REAL(wp) ::  e_int_l            !:
    121128    REAL(wp) ::  e_int_u            !:
     
    123130    REAL(wp) ::  exp_arg            !:
    124131    REAL(wp) ::  exp_term           !:
    125     REAL(wp) ::  fs_int             !:
    126132    REAL(wp) ::  gg                 !:
    127133    REAL(wp) ::  height_int         !:
     
    129135    REAL(wp) ::  lagr_timescale     !:
    130136    REAL(wp) ::  location(1:30,1:3) !:
    131     REAL(wp) ::  log_z_z0_int       !:
    132137    REAL(wp) ::  random_gauss       !:
    133     REAL(wp) ::  u_int              !:
    134138    REAL(wp) ::  u_int_l            !:
    135139    REAL(wp) ::  u_int_u            !:
    136140    REAL(wp) ::  us_int             !:
    137     REAL(wp) ::  v_int              !:
    138141    REAL(wp) ::  v_int_l            !:
    139142    REAL(wp) ::  v_int_u            !:
    140143    REAL(wp) ::  vv_int             !:
    141     REAL(wp) ::  w_int              !:
    142144    REAL(wp) ::  w_int_l            !:
    143145    REAL(wp) ::  w_int_u            !:
     
    153155    REAL(wp), DIMENSION(1:30) ::  ei      !:
    154156
     157    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !:
     158    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !:
     159    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !:
     160    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !:
     161    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !:
     162    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !:
     163    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !:
     164    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !:
     165    REAL(wp), DIMENSION(number_of_particles) ::  log_z_z0_int !:
     166    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !:
     167    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !:
     168    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !:
     169    REAL(wp), DIMENSION(number_of_particles) ::  xv           !:
     170    REAL(wp), DIMENSION(number_of_particles) ::  yv           !:
     171    REAL(wp), DIMENSION(number_of_particles) ::  zv           !:
     172
     173    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !:
     174
     175    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
     176
    155177!
    156178!-- Determine height of Prandtl layer and distance between Prandtl-layer
     
    159181!-- (for particles below first vertical grid level).
    160182    z_p      = zu(nzb+1) - zw(nzb)
    161     d_z_p_z0 = 1.0 / ( z_p - z0_av_global )
    162 
    163     DO  n = 1, number_of_particles
    164 
    165 !
    166 !--    Move particle only if the LES timestep has not (approximately) been
    167 !--    reached
    168        IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 )  CYCLE
    169 !
    170 !--    Determine bottom index
    171        k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    172               + offset_ocean_nzt       ! only exact if equidistant
    173 !
    174 !--    Interpolation of the u velocity component onto particle position. 
    175 !--    Particles are interpolation bi-linearly in the horizontal and a
    176 !--    linearly in the vertical. An exception is made for particles below
    177 !--    the first vertical grid level in case of a prandtl layer. In this
    178 !--    case the horizontal particle velocity components are determined using
    179 !--    Monin-Obukhov relations (if branch).
    180 !--    First, check if particle is located below first vertical grid level
    181 !--    (Prandtl-layer height)
    182        IF ( prandtl_layer .AND. particles(n)%z < z_p )  THEN
    183 !
    184 !--       Resolved-scale horizontal particle velocity is zero below z0.
    185           IF ( particles(n)%z < z0_av_global )  THEN
    186 
    187              u_int = 0.0
    188 
     183    d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
     184
     185    start_index = grid_particles(kp,jp,ip)%start_index
     186    end_index   = grid_particles(kp,jp,ip)%end_index
     187
     188    xv = particles(1:number_of_particles)%x
     189    yv = particles(1:number_of_particles)%y
     190    zv = particles(1:number_of_particles)%z
     191
     192    DO  nb = 0, 7
     193
     194       i = ip
     195       j = jp + block_offset(nb)%j_off
     196       k = kp + block_offset(nb)%k_off
     197
     198!
     199!--    Interpolate u velocity-component
     200       DO  n = start_index(nb), end_index(nb)
     201!
     202!--       Interpolation of the u velocity component onto particle position. 
     203!--       Particles are interpolation bi-linearly in the horizontal and a
     204!--       linearly in the vertical. An exception is made for particles below
     205!--       the first vertical grid level in case of a prandtl layer. In this
     206!--       case the horizontal particle velocity components are determined using
     207!--       Monin-Obukhov relations (if branch).
     208!--       First, check if particle is located below first vertical grid level
     209!--       (Prandtl-layer height)
     210          IF ( prandtl_layer  .AND.  particles(n)%z < z_p )  THEN
     211!
     212!--          Resolved-scale horizontal particle velocity is zero below z0.
     213             IF ( particles(n)%z < z0_av_global )  THEN
     214                u_int(n) = 0.0_wp
     215             ELSE
     216!
     217!--             Determine the sublayer. Further used as index.
     218                height_p     = ( particles(n)%z - z0_av_global )            &
     219                                     * REAL( number_of_sublayers, KIND=wp ) &
     220                                     * d_z_p_z0
     221!
     222!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
     223!--             interpolate linearly between precalculated logarithm.
     224                log_z_z0_int(n) = log_z_z0(INT(height_p))                      &
     225                                 + ( height_p - INT(height_p) )                &
     226                                 * ( log_z_z0(INT(height_p)+1)                 &
     227                                      - log_z_z0(INT(height_p))                &
     228                                   )
     229!
     230!--             Neutral solution is applied for all situations, e.g. also for
     231!--             unstable and stable situations. Even though this is not exact
     232!--             this saves a lot of CPU time since several calls of intrinsic
     233!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     234!--             as sensitivity studies revealed no significant effect of
     235!--             using the neutral solution also for un/stable situations.
     236!--             Calculated left and bottom index on u grid.
     237                us_int = 0.5_wp * ( us(j,i) + us(j,i-1) ) 
     238
     239                u_int  = -usws(j,i) / ( us_int * kappa + 1E-10_wp )           &
     240                         * log_z_z0_int(n)
     241
     242             ENDIF
     243!
     244!--       Particle above the first grid level. Bi-linear interpolation in the
     245!--       horizontal and linear interpolation in the vertical direction.
    189246          ELSE
    190 !
    191 !--          Determine the sublayer. Further used as index.
    192              height_p     = ( particles(n)%z - z0_av_global )            &
    193                                   * REAL( number_of_sublayers, KIND=wp ) &
    194                                   * d_z_p_z0
    195 !
    196 !--          Calculate LOG(z/z0) for exact particle height. Therefore,   
    197 !--          interpolate linearly between precalculated logarithm.
    198              log_z_z0_int = log_z_z0(INT(height_p))                     &
    199                           + ( height_p - INT(height_p) )                &
    200                           * ( log_z_z0(INT(height_p)+1)                 &
    201                                - log_z_z0(INT(height_p))                &
    202                             )
    203 !
    204 !--          Neutral solution is applied for all situations, e.g. also for
    205 !--          unstable and stable situations. Even though this is not exact
    206 !--          this saves a lot of CPU time since several calls of intrinsic
    207 !--          FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    208 !--          as sensitivity studies revealed no significant effect of
    209 !--          using the neutral solution also for un/stable situations.
    210 !--          Calculated left and bottom index on u grid.
    211              i      = ( particles(n)%x + 0.5 * dx ) * ddx
    212              j      =   particles(n)%y              * ddy
    213 
    214              us_int = 0.5 * ( us(j,i) + us(j,i-1) ) 
    215 
    216              u_int  = -usws(j,i) / ( us_int * kappa + 1E-10_wp )           &
    217                       * log_z_z0_int
    218 
    219           ENDIF
    220 !
    221 !--    Particle above the first grid level. Bi-linear interpolation in the
    222 !--    horizontal and linear interpolation in the vertical direction.
    223        ELSE
    224 !
    225 !--       Interpolate u velocity-component, determine left, front, bottom
    226 !--       index of u-array. Adopt k index from above
    227           i = ( particles(n)%x + 0.5 * dx ) * ddx
    228           j =   particles(n)%y * ddy
    229 !
    230 !--       Interpolation of the velocity components in the xy-plane
    231           x  = particles(n)%x + ( 0.5 - i ) * dx
    232           y  = particles(n)%y - j * dy
    233           aa = x**2          + y**2
    234           bb = ( dx - x )**2 + y**2
    235           cc = x**2          + ( dy - y )**2
    236           dd = ( dx - x )**2 + ( dy - y )**2
    237           gg = aa + bb + cc + dd
    238 
    239           u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)  &
    240                     + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1)&
    241                     ) / ( 3.0 * gg ) - u_gtrans
    242 
    243           IF ( k+1 == nzt+1 )  THEN
    244 
    245              u_int = u_int_l
    246 
    247           ELSE
    248 
    249              u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
    250                        + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
    251                        ) / ( 3.0 * gg ) - u_gtrans
    252 
    253              u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz *            &
    254                             ( u_int_u - u_int_l )
    255 
    256           ENDIF
    257 
    258        ENDIF
    259 
    260 !
    261 !--    Same procedure for interpolation of the v velocity-component.
    262        IF ( prandtl_layer .AND. particles(n)%z < z_p )  THEN
    263 !
    264 !--       Resolved-scale horizontal particle velocity is zero below z0.
    265           IF ( particles(n)%z < z0_av_global )  THEN
    266 
    267              v_int = 0.0
    268 
    269           ELSE       
    270 !
    271 !--          Neutral solution is applied for all situations, e.g. also for
    272 !--          unstable and stable situations. Even though this is not exact
    273 !--          this saves a lot of CPU time since several calls of intrinsic
    274 !--          FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    275 !--          as sensitivity studies revealed no significant effect of
    276 !--          using the neutral solution also for un/stable situations.
    277 !--          Calculated left and bottom index on v grid.
    278               i      =   particles(n)%x              * ddx
    279               j      = ( particles(n)%y + 0.5 * dy ) * ddy
    280 
    281               us_int = 0.5 * ( us(j,i) + us(j-1,i) )   
    282 
    283               v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10_wp )             &
    284                        * log_z_z0_int
    285 
    286           ENDIF
    287 !
    288 !--    Particle above the first grid level. Bi-linear interpolation in the
    289 !--    horizontal and linear interpolation in the vertical direction.
    290        ELSE
    291           i =   particles(n)%x * ddx
    292           j = ( particles(n)%y + 0.5 * dy ) * ddy
    293           x  = particles(n)%x - i * dx
    294           y  = particles(n)%y + ( 0.5 - j ) * dy
    295           aa = x**2          + y**2
    296           bb = ( dx - x )**2 + y**2
    297           cc = x**2          + ( dy - y )**2
    298           dd = ( dx - x )**2 + ( dy - y )**2
    299           gg = aa + bb + cc + dd
    300 
    301           v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)  &
    302                     + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1)&
    303                     ) / ( 3.0 * gg ) - v_gtrans
    304           IF ( k+1 == nzt+1 )  THEN
    305              v_int = v_int_l
    306           ELSE
    307              v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
    308                        + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
    309                        ) / ( 3.0 * gg ) - v_gtrans
    310              v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz *           &
    311                             ( v_int_u - v_int_l )
    312           ENDIF
    313 
    314        ENDIF
    315 
    316 !
    317 !--    Same procedure for interpolation of the w velocity-component
    318        IF ( vertical_particle_advection(particles(n)%group) )  THEN
    319           i = particles(n)%x * ddx
    320           j = particles(n)%y * ddy
    321           k = particles(n)%z / dz + offset_ocean_nzt_m1
    322 
    323           x  = particles(n)%x - i * dx
    324           y  = particles(n)%y - j * dy
    325           aa = x**2          + y**2
    326           bb = ( dx - x )**2 + y**2
    327           cc = x**2          + ( dy - y )**2
    328           dd = ( dx - x )**2 + ( dy - y )**2
    329           gg = aa + bb + cc + dd
    330 
    331           w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)    &
    332                     + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1)  &
    333                        ) / ( 3.0 * gg )
    334           IF ( k+1 == nzt+1 )  THEN
    335              w_int = w_int_l
    336           ELSE
    337              w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
    338                          ( gg-bb ) * w(k+1,j,i+1) + &
    339                          ( gg-cc ) * w(k+1,j+1,i) + &
    340                          ( gg-dd ) * w(k+1,j+1,i+1) &
    341                         ) / ( 3.0 * gg )
    342              w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
    343                                ( w_int_u - w_int_l )
    344           ENDIF
    345        ELSE
    346           w_int = 0.0
    347        ENDIF
    348 
    349 !
    350 !--    Interpolate and calculate quantities needed for calculating the SGS
    351 !--    velocities
    352        IF ( use_sgs_for_particles )  THEN
    353 !
    354 !--       Interpolate TKE
    355           i = particles(n)%x * ddx
    356           j = particles(n)%y * ddy
    357           k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    358               + offset_ocean_nzt                      ! only exact if eq.dist
    359 
    360           IF ( topography == 'flat' )  THEN
    361 
    362              x  = particles(n)%x - i * dx
    363              y  = particles(n)%y - j * dy
     247
     248             x  = xv(n) + ( 0.5_wp - i ) * dx
     249             y  = yv(n) - j * dy
    364250             aa = x**2          + y**2
    365251             bb = ( dx - x )**2 + y**2
     
    368254             gg = aa + bb + cc + dd
    369255
    370              e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
    371                        + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
    372                        ) / ( 3.0 * gg )
    373 
    374              IF ( k+1 == nzt+1 )  THEN
    375                 e_int = e_int_l
     256             u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
     257                         + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
     258                         u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
     259
     260             IF ( k == nzt )  THEN
     261                u_int(n) = u_int_l
    376262             ELSE
    377                 e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
    378                             ( gg - bb ) * e(k+1,j,i+1) + &
    379                             ( gg - cc ) * e(k+1,j+1,i) + &
    380                             ( gg - dd ) * e(k+1,j+1,i+1) &
    381                           ) / ( 3.0 * gg )
    382                 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
    383                                   ( e_int_u - e_int_l )
    384              ENDIF
    385 
    386 !
    387 !--          Interpolate the TKE gradient along x (adopt incides i,j,k and
    388 !--          all position variables from above (TKE))
    389              de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
    390                              ( gg - bb ) * de_dx(k,j,i+1) + &
    391                              ( gg - cc ) * de_dx(k,j+1,i) + &
    392                              ( gg - dd ) * de_dx(k,j+1,i+1) &
    393                            ) / ( 3.0 * gg )
    394 
    395              IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    396                 de_dx_int = de_dx_int_l
     263                u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
     264                            + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
     265                            u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
     266                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *                  &
     267                           ( u_int_u - u_int_l )
     268             ENDIF
     269          ENDIF
     270
     271       ENDDO
     272
     273       i = ip + block_offset(nb)%i_off
     274       j = jp
     275       k = kp + block_offset(nb)%k_off
     276!
     277!--    Same procedure for interpolation of the v velocity-component
     278       DO  n = start_index(nb), end_index(nb)
     279          IF ( prandtl_layer  .AND.  particles(n)%z < z_p )  THEN
     280
     281             IF ( particles(n)%z < z0_av_global )  THEN
     282!
     283!--             Resolved-scale horizontal particle velocity is zero below z0.
     284                v_int(n) = 0.0_wp
     285             ELSE       
     286!
     287!--             Neutral solution is applied for all situations, e.g. also for
     288!--             unstable and stable situations. Even though this is not exact
     289!--             this saves a lot of CPU time since several calls of intrinsic
     290!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     291!--             as sensitivity studies revealed no significant effect of
     292!--             using the neutral solution also for un/stable situations.
     293!--             Calculated left and bottom index on v grid.
     294                us_int = 0.5_wp * ( us(j,i) + us(j-1,i) )   
     295
     296                v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10_wp )             &
     297                         * log_z_z0_int(n)
     298             ENDIF
     299          ELSE
     300             x  = xv(n) - i * dx
     301             y  = yv(n) + ( 0.5_wp - j ) * dy
     302             aa = x**2          + y**2
     303             bb = ( dx - x )**2 + y**2
     304             cc = x**2          + ( dy - y )**2
     305             dd = ( dx - x )**2 + ( dy - y )**2
     306             gg = aa + bb + cc + dd
     307
     308             v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
     309                       + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
     310                       ) / ( 3.0_wp * gg ) - v_gtrans
     311
     312             IF ( k == nzt )  THEN
     313                v_int(n) = v_int_l
    397314             ELSE
    398                 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
    399                                 ( gg - bb ) * de_dx(k+1,j,i+1) + &
    400                                 ( gg - cc ) * de_dx(k+1,j+1,i) + &
    401                                 ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    402                               ) / ( 3.0 * gg )
    403                 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
    404                                           ( de_dx_int_u - de_dx_int_l )
    405              ENDIF
    406 
    407 !
    408 !--          Interpolate the TKE gradient along y
    409              de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
    410                              ( gg - bb ) * de_dy(k,j,i+1) + &
    411                              ( gg - cc ) * de_dy(k,j+1,i) + &
    412                              ( gg - dd ) * de_dy(k,j+1,i+1) &
    413                            ) / ( 3.0 * gg )
    414              IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    415                 de_dy_int = de_dy_int_l
     315                v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
     316                          + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
     317                          ) / ( 3.0_wp * gg ) - v_gtrans
     318                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
     319                                  ( v_int_u - v_int_l )
     320             ENDIF
     321          ENDIF
     322
     323       ENDDO
     324
     325       i = ip + block_offset(nb)%i_off
     326       j = jp + block_offset(nb)%j_off
     327       k = kp-1
     328!
     329!--    Same procedure for interpolation of the w velocity-component
     330       DO  n = start_index(nb), end_index(nb)
     331
     332          IF ( vertical_particle_advection(particles(n)%group) )  THEN
     333
     334             x  = xv(n) - i * dx
     335             y  = yv(n) - j * dy
     336             aa = x**2          + y**2
     337             bb = ( dx - x )**2 + y**2
     338             cc = x**2          + ( dy - y )**2
     339             dd = ( dx - x )**2 + ( dy - y )**2
     340             gg = aa + bb + cc + dd
     341
     342             w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
     343                       + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
     344                       ) / ( 3.0_wp * gg )
     345
     346             IF ( k == nzt )  THEN
     347                w_int(n) = w_int_l
    416348             ELSE
    417                 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    418                                 ( gg - bb ) * de_dy(k+1,j,i+1) + &
    419                                 ( gg - cc ) * de_dy(k+1,j+1,i) + &
    420                                 ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    421                               ) / ( 3.0 * gg )
    422                 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
    423                                           ( de_dy_int_u - de_dy_int_l )
    424              ENDIF
    425 
    426 !
    427 !--          Interpolate the TKE gradient along z
    428              IF ( particles(n)%z < 0.5 * dz ) THEN
    429                 de_dz_int = 0.0
    430              ELSE
    431                 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
    432                                 ( gg - bb ) * de_dz(k,j,i+1) + &
    433                                 ( gg - cc ) * de_dz(k,j+1,i) + &
    434                                 ( gg - dd ) * de_dz(k,j+1,i+1) &
    435                               ) / ( 3.0 * gg )
     349                w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
     350                            ( gg-bb ) * w(k+1,j,i+1) + &
     351                            ( gg-cc ) * w(k+1,j+1,i) + &
     352                            ( gg-dd ) * w(k+1,j+1,i+1) &
     353                          ) / ( 3.0_wp * gg )
     354                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
     355                           ( w_int_u - w_int_l )
     356             ENDIF
     357
     358          ELSE
     359
     360             w_int(n) = 0.0_wp
     361
     362          ENDIF
     363
     364       ENDDO
     365
     366    ENDDO
     367
     368!-- Interpolate and calculate quantities needed for calculating the SGS
     369!-- velocities
     370    IF ( use_sgs_for_particles )  THEN
     371
     372       IF ( topography == 'flat' )  THEN
     373
     374          DO  nb = 0,7
     375
     376             i = ip + block_offset(nb)%i_off
     377             j = jp + block_offset(nb)%j_off
     378             k = kp + block_offset(nb)%k_off
     379
     380             DO  n = start_index(nb), end_index(nb)
     381!
     382!--             Interpolate TKE
     383                x  = xv(n) - i * dx
     384                y  = yv(n) - j * dy
     385                aa = x**2          + y**2
     386                bb = ( dx - x )**2 + y**2
     387                cc = x**2          + ( dy - y )**2
     388                dd = ( dx - x )**2 + ( dy - y )**2
     389                gg = aa + bb + cc + dd
     390
     391                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
     392                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
     393                          ) / ( 3.0_wp * gg )
     394
     395                IF ( k+1 == nzt+1 )  THEN
     396                   e_int(n) = e_int_l
     397                ELSE
     398                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
     399                               ( gg - bb ) * e(k+1,j,i+1) + &
     400                               ( gg - cc ) * e(k+1,j+1,i) + &
     401                               ( gg - dd ) * e(k+1,j+1,i+1) &
     402                            ) / ( 3.0_wp * gg )
     403                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
     404                                     ( e_int_u - e_int_l )
     405                ENDIF
     406!
     407!--             Needed to avoid NaN particle velocities
     408                IF ( e_int(n) == 0.0_wp )  THEN
     409                   e_int(n) = 1.0E-20_wp
     410                ENDIF
     411!
     412!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
     413!--             all position variables from above (TKE))
     414                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
     415                                ( gg - bb ) * de_dx(k,j,i+1) + &
     416                                ( gg - cc ) * de_dx(k,j+1,i) + &
     417                                ( gg - dd ) * de_dx(k,j+1,i+1) &
     418                               ) / ( 3.0_wp * gg )
    436419
    437420                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    438                    de_dz_int = de_dz_int_l
     421                   de_dx_int(n) = de_dx_int_l
    439422                ELSE
    440                    de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
    441                                    ( gg - bb ) * de_dz(k+1,j,i+1) + &
    442                                    ( gg - cc ) * de_dz(k+1,j+1,i) + &
    443                                    ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    444                                  ) / ( 3.0 * gg )
    445                    de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
    446                                              ( de_dz_int_u - de_dz_int_l )
    447                 ENDIF
    448              ENDIF
    449 
    450 !
    451 !--          Interpolate the dissipation of TKE
    452              diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
    453                             ( gg - bb ) * diss(k,j,i+1) + &
    454                             ( gg - cc ) * diss(k,j+1,i) + &
    455                             ( gg - dd ) * diss(k,j+1,i+1) &
    456                            ) / ( 3.0 * gg )
    457 
    458              IF ( k+1 == nzt+1 )  THEN
    459                 diss_int = diss_int_l
    460              ELSE
    461                 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
    462                                ( gg - bb ) * diss(k+1,j,i+1) + &
    463                                ( gg - cc ) * diss(k+1,j+1,i) + &
    464                                ( gg - dd ) * diss(k+1,j+1,i+1) &
    465                               ) / ( 3.0 * gg )
    466                 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
    467                                         ( diss_int_u - diss_int_l )
    468              ENDIF
    469 
    470           ELSE
    471 
     423                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
     424                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
     425                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
     426                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
     427                                  ) / ( 3.0_wp * gg )
     428                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
     429                                              ( de_dx_int_u - de_dx_int_l )
     430                ENDIF
     431!
     432!--             Interpolate the TKE gradient along y
     433                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
     434                                ( gg - bb ) * de_dy(k,j,i+1) + &
     435                                ( gg - cc ) * de_dy(k,j+1,i) + &
     436                                ( gg - dd ) * de_dy(k,j+1,i+1) &
     437                               ) / ( 3.0_wp * gg )
     438                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     439                   de_dy_int(n) = de_dy_int_l
     440                ELSE
     441                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
     442                                  ( gg - bb ) * de_dy(k+1,j,i+1) + &
     443                                  ( gg - cc ) * de_dy(k+1,j+1,i) + &
     444                                  ( gg - dd ) * de_dy(k+1,j+1,i+1) &
     445                                  ) / ( 3.0_wp * gg )
     446                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
     447                                              ( de_dy_int_u - de_dy_int_l )
     448                ENDIF
     449
     450!
     451!--             Interpolate the TKE gradient along z
     452                IF ( zv(n) < 0.5_wp * dz )  THEN
     453                   de_dz_int(n) = 0.0_wp
     454                ELSE
     455                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
     456                                   ( gg - bb ) * de_dz(k,j,i+1) + &
     457                                   ( gg - cc ) * de_dz(k,j+1,i) + &
     458                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
     459                                  ) / ( 3.0_wp * gg )
     460
     461                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     462                      de_dz_int(n) = de_dz_int_l
     463                   ELSE
     464                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
     465                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
     466                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
     467                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
     468                                     ) / ( 3.0_wp * gg )
     469                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
     470                                                 ( de_dz_int_u - de_dz_int_l )
     471                   ENDIF
     472                ENDIF
     473
     474!
     475!--             Interpolate the dissipation of TKE
     476                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
     477                               ( gg - bb ) * diss(k,j,i+1) + &
     478                               ( gg - cc ) * diss(k,j+1,i) + &
     479                               ( gg - dd ) * diss(k,j+1,i+1) &
     480                              ) / ( 3.0_wp * gg )
     481
     482                IF ( k == nzt )  THEN
     483                   diss_int(n) = diss_int_l
     484                ELSE
     485                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
     486                                  ( gg - bb ) * diss(k+1,j,i+1) + &
     487                                  ( gg - cc ) * diss(k+1,j+1,i) + &
     488                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
     489                                 ) / ( 3.0_wp * gg )
     490                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
     491                                           ( diss_int_u - diss_int_l )
     492                ENDIF
     493
     494             ENDDO
     495          ENDDO
     496
     497       ELSE ! non-flat topography, e.g., buildings
     498
     499          DO  n = 1, number_of_particles
     500
     501             i = particles(n)%x * ddx
     502             j = particles(n)%y * ddy
     503             k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
     504                 + offset_ocean_nzt                      ! only exact if eq.dist
    472505!
    473506!--          In case that there are buildings it has to be determined
     
    484517
    485518             gp_outside_of_building = 0
    486              location = 0.0
     519             location = 0.0_wp
    487520             num_gp = 0
    488521
     
    492525                location(num_gp,1) = i * dx
    493526                location(num_gp,2) = j * dy
    494                 location(num_gp,3) = k * dz - 0.5 * dz
     527                location(num_gp,3) = k * dz - 0.5_wp * dz
    495528                ei(num_gp)     = e(k,j,i)
    496529                dissi(num_gp)  = diss(k,j,i)
     
    506539                location(num_gp,1) = i * dx
    507540                location(num_gp,2) = (j+1) * dy
    508                 location(num_gp,3) = k * dz - 0.5 * dz
     541                location(num_gp,3) = k * dz - 0.5_wp * dz
    509542                ei(num_gp)     = e(k,j+1,i)
    510543                dissi(num_gp)  = diss(k,j+1,i)
     
    519552                location(num_gp,1) = i * dx
    520553                location(num_gp,2) = j * dy
    521                 location(num_gp,3) = (k+1) * dz - 0.5 * dz
     554                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    522555                ei(num_gp)     = e(k+1,j,i)
    523556                dissi(num_gp)  = diss(k+1,j,i)
     
    533566                location(num_gp,1) = i * dx
    534567                location(num_gp,2) = (j+1) * dy
    535                 location(num_gp,3) = (k+1) * dz - 0.5 * dz
     568                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    536569                ei(num_gp)     = e(k+1,j+1,i)
    537570                dissi(num_gp)  = diss(k+1,j+1,i)
     
    547580                location(num_gp,1) = (i+1) * dx
    548581                location(num_gp,2) = j * dy
    549                 location(num_gp,3) = k * dz - 0.5 * dz
     582                location(num_gp,3) = k * dz - 0.5_wp * dz
    550583                ei(num_gp)     = e(k,j,i+1)
    551584                dissi(num_gp)  = diss(k,j,i+1)
     
    555588             ENDIF
    556589
    557              IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
     590             IF ( k > nzb_s_inner(j+1,i+1)  .OR. nzb_s_inner(j+1,i+1) == 0 ) &
    558591             THEN
    559592                num_gp = num_gp + 1
     
    561594                location(num_gp,1) = (i+1) * dx
    562595                location(num_gp,2) = (j+1) * dy
    563                 location(num_gp,3) = k * dz - 0.5 * dz
     596                location(num_gp,3) = k * dz - 0.5_wp * dz
    564597                ei(num_gp)     = e(k,j+1,i+1)
    565598                dissi(num_gp)  = diss(k,j+1,i+1)
     
    575608                location(num_gp,1) = (i+1) * dx
    576609                location(num_gp,2) = j * dy
    577                 location(num_gp,3) = (k+1) * dz - 0.5 * dz
     610                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    578611                ei(num_gp)     = e(k+1,j,i+1)
    579612                dissi(num_gp)  = diss(k+1,j,i+1)
     
    583616             ENDIF
    584617
    585              IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
     618             IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR. nzb_s_inner(j+1,i+1) == 0)&
    586619             THEN
    587620                num_gp = num_gp + 1
     
    589622                location(num_gp,1) = (i+1) * dx
    590623                location(num_gp,2) = (j+1) * dy
    591                 location(num_gp,3) = (k+1) * dz - 0.5 * dz
     624                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    592625                ei(num_gp)     = e(k+1,j+1,i+1)
    593626                dissi(num_gp)  = diss(k+1,j+1,i+1)
     
    610643                gg = aa + bb + cc + dd
    611644
    612                 e_int_l = (( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1) &
    613                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
    614                           ) / ( 3.0 * gg )
    615 
    616                 IF ( k+1 == nzt+1 )  THEN
    617                    e_int = e_int_l
     645                e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)  &
     646                          + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
     647                          ) / ( 3.0_wp * gg )
     648
     649                IF ( k == nzt )  THEN
     650                   e_int(n) = e_int_l
    618651                ELSE
    619652                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
     
    621654                               ( gg - cc ) * e(k+1,j+1,i) + &
    622655                               ( gg - dd ) * e(k+1,j+1,i+1) &
    623                              ) / ( 3.0 * gg )
    624                    e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
     656                             ) / ( 3.0_wp * gg )
     657                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
    625658                                     ( e_int_u - e_int_l )
    626659                ENDIF
    627 
     660!
     661!--             Needed to avoid NaN particle velocities
     662                IF ( e_int(n) == 0.0_wp )  THEN
     663                   e_int(n) = 1.0E-20_wp
     664                ENDIF
    628665!
    629666!--             Interpolate the TKE gradient along x (adopt incides i,j,k
     
    633670                                ( gg - cc ) * de_dx(k,j+1,i) + &
    634671                                ( gg - dd ) * de_dx(k,j+1,i+1) &
    635                               ) / ( 3.0 * gg )
    636 
    637                 IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    638                    de_dx_int = de_dx_int_l
     672                              ) / ( 3.0_wp * gg )
     673
     674                IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
     675                   de_dx_int(n) = de_dx_int_l
    639676                ELSE
    640677                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
     
    642679                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
    643680                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    644                                  ) / ( 3.0 * gg )
    645                    de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
     681                                 ) / ( 3.0_wp * gg )
     682                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
    646683                                           dz * ( de_dx_int_u - de_dx_int_l )
    647684                ENDIF
     
    653690                                ( gg - cc ) * de_dy(k,j+1,i) + &
    654691                                ( gg - dd ) * de_dy(k,j+1,i+1) &
    655                               ) / ( 3.0 * gg )
     692                              ) / ( 3.0_wp * gg )
    656693                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    657                    de_dy_int = de_dy_int_l
     694                   de_dy_int(n) = de_dy_int_l
    658695                ELSE
    659696                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
     
    661698                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
    662699                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    663                                  ) / ( 3.0 * gg )
    664                    de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
     700                                 ) / ( 3.0_wp * gg )
     701                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
    665702                                           dz * ( de_dy_int_u - de_dy_int_l )
    666703                ENDIF
     
    668705!
    669706!--             Interpolate the TKE gradient along z
    670                 IF ( particles(n)%z < 0.5 * dz )  THEN
    671                    de_dz_int = 0.0
     707                IF ( zv(n) < 0.5_wp * dz )  THEN
     708                   de_dz_int(n) = 0.0_wp
    672709                ELSE
    673710                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
     
    675712                                   ( gg - cc ) * de_dz(k,j+1,i) + &
    676713                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
    677                                  ) / ( 3.0 * gg )
     714                                 ) / ( 3.0_wp * gg )
    678715
    679716                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    680                       de_dz_int = de_dz_int_l
     717                      de_dz_int(n) = de_dz_int_l
    681718                   ELSE
    682719                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
     
    684721                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
    685722                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    686                                     ) / ( 3.0 * gg )
    687                       de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
     723                                    ) / ( 3.0_wp * gg )
     724                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
    688725                                           dz * ( de_dz_int_u - de_dz_int_l )
    689726                   ENDIF
     
    696733                               ( gg - cc ) * diss(k,j+1,i) + &
    697734                               ( gg - dd ) * diss(k,j+1,i+1) &
    698                              ) / ( 3.0 * gg )
    699 
    700                 IF ( k+1 == nzt+1 )  THEN
    701                    diss_int = diss_int_l
     735                             ) / ( 3.0_wp * gg )
     736
     737                IF ( k == nzt )  THEN
     738                   diss_int(n) = diss_int_l
    702739                ELSE
    703740                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
     
    705742                                  ( gg - cc ) * diss(k+1,j+1,i) + &
    706743                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
    707                                 ) / ( 3.0 * gg )
    708                    diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
     744                                ) / ( 3.0_wp * gg )
     745                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
    709746                                           ( diss_int_u - diss_int_l )
    710747                ENDIF
     
    718755                     gp_outside_of_building(5) == 0 )  THEN
    719756                   num_gp = num_gp + 1
    720                    location(num_gp,1) = i * dx + 0.5 * dx
     757                   location(num_gp,1) = i * dx + 0.5_wp * dx
    721758                   location(num_gp,2) = j * dy
    722                    location(num_gp,3) = k * dz - 0.5 * dz
     759                   location(num_gp,3) = k * dz - 0.5_wp * dz
    723760                   ei(num_gp)     = e(k,j,i)
    724761                   dissi(num_gp)  = diss(k,j,i)
    725                    de_dxi(num_gp) = 0.0
     762                   de_dxi(num_gp) = 0.0_wp
    726763                   de_dyi(num_gp) = de_dy(k,j,i)
    727764                   de_dzi(num_gp) = de_dz(k,j,i)
     
    731768                   gp_outside_of_building(1) == 0 )  THEN
    732769                   num_gp = num_gp + 1
    733                    location(num_gp,1) = i * dx + 0.5 * dx
     770                   location(num_gp,1) = i * dx + 0.5_wp * dx
    734771                   location(num_gp,2) = j * dy
    735                    location(num_gp,3) = k * dz - 0.5 * dz
     772                   location(num_gp,3) = k * dz - 0.5_wp * dz
    736773                   ei(num_gp)     = e(k,j,i+1)
    737774                   dissi(num_gp)  = diss(k,j,i+1)
    738                    de_dxi(num_gp) = 0.0
     775                   de_dxi(num_gp) = 0.0_wp
    739776                   de_dyi(num_gp) = de_dy(k,j,i+1)
    740777                   de_dzi(num_gp) = de_dz(k,j,i+1)
     
    748785                   num_gp = num_gp + 1
    749786                   location(num_gp,1) = (i+1) * dx
    750                    location(num_gp,2) = j * dy + 0.5 * dy
    751                    location(num_gp,3) = k * dz - 0.5 * dz
     787                   location(num_gp,2) = j * dy + 0.5_wp * dy
     788                   location(num_gp,3) = k * dz - 0.5_wp * dz
    752789                   ei(num_gp)     = e(k,j,i+1)
    753790                   dissi(num_gp)  = diss(k,j,i+1)
    754791                   de_dxi(num_gp) = de_dx(k,j,i+1)
    755                    de_dyi(num_gp) = 0.0
     792                   de_dyi(num_gp) = 0.0_wp
    756793                   de_dzi(num_gp) = de_dz(k,j,i+1)
    757794                ENDIF
     
    761798                   num_gp = num_gp + 1
    762799                   location(num_gp,1) = (i+1) * dx
    763                    location(num_gp,2) = j * dy + 0.5 * dy
    764                    location(num_gp,3) = k * dz - 0.5 * dz
     800                   location(num_gp,2) = j * dy + 0.5_wp * dy
     801                   location(num_gp,3) = k * dz - 0.5_wp * dz
    765802                   ei(num_gp)     = e(k,j+1,i+1)
    766803                   dissi(num_gp)  = diss(k,j+1,i+1)
    767804                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
    768                    de_dyi(num_gp) = 0.0
     805                   de_dyi(num_gp) = 0.0_wp
    769806                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
    770807                ENDIF
     
    776813                     gp_outside_of_building(6) == 0 )  THEN
    777814                   num_gp = num_gp + 1
    778                    location(num_gp,1) = i * dx + 0.5 * dx
     815                   location(num_gp,1) = i * dx + 0.5_wp * dx
    779816                   location(num_gp,2) = (j+1) * dy
    780                    location(num_gp,3) = k * dz - 0.5 * dz
     817                   location(num_gp,3) = k * dz - 0.5_wp * dz
    781818                   ei(num_gp)     = e(k,j+1,i)
    782819                   dissi(num_gp)  = diss(k,j+1,i)
    783                    de_dxi(num_gp) = 0.0
     820                   de_dxi(num_gp) = 0.0_wp
    784821                   de_dyi(num_gp) = de_dy(k,j+1,i)
    785822                   de_dzi(num_gp) = de_dz(k,j+1,i)
     
    789826                     gp_outside_of_building(2) == 0 )  THEN
    790827                   num_gp = num_gp + 1
    791                    location(num_gp,1) = i * dx + 0.5 * dx
     828                   location(num_gp,1) = i * dx + 0.5_wp * dx
    792829                   location(num_gp,2) = (j+1) * dy
    793                    location(num_gp,3) = k * dz - 0.5 * dz
     830                   location(num_gp,3) = k * dz - 0.5_wp * dz
    794831                   ei(num_gp)     = e(k,j+1,i+1)
    795832                   dissi(num_gp)  = diss(k,j+1,i+1)
    796                    de_dxi(num_gp) = 0.0
     833                   de_dxi(num_gp) = 0.0_wp
    797834                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
    798835                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
     
    806843                   num_gp = num_gp + 1
    807844                   location(num_gp,1) = i * dx
    808                    location(num_gp,2) = j * dy + 0.5 * dy
    809                    location(num_gp,3) = k * dz - 0.5 * dz
     845                   location(num_gp,2) = j * dy + 0.5_wp * dy
     846                   location(num_gp,3) = k * dz - 0.5_wp * dz
    810847                   ei(num_gp)     = e(k,j,i)
    811848                   dissi(num_gp)  = diss(k,j,i)
    812849                   de_dxi(num_gp) = de_dx(k,j,i)
    813                    de_dyi(num_gp) = 0.0
     850                   de_dyi(num_gp) = 0.0_wp
    814851                   de_dzi(num_gp) = de_dz(k,j,i)
    815852                ENDIF
     
    819856                   num_gp = num_gp + 1
    820857                   location(num_gp,1) = i * dx
    821                    location(num_gp,2) = j * dy + 0.5 * dy
    822                    location(num_gp,3) = k * dz - 0.5 * dz
     858                   location(num_gp,2) = j * dy + 0.5_wp * dy
     859                   location(num_gp,3) = k * dz - 0.5_wp * dz
    823860                   ei(num_gp)     = e(k,j+1,i)
    824861                   dissi(num_gp)  = diss(k,j+1,i)
    825862                   de_dxi(num_gp) = de_dx(k,j+1,i)
    826                    de_dyi(num_gp) = 0.0
     863                   de_dyi(num_gp) = 0.0_wp
    827864                   de_dzi(num_gp) = de_dz(k,j+1,i)
    828865                ENDIF
     
    834871                     gp_outside_of_building(7) == 0 )  THEN
    835872                   num_gp = num_gp + 1
    836                    location(num_gp,1) = i * dx + 0.5 * dx
     873                   location(num_gp,1) = i * dx + 0.5_wp * dx
    837874                   location(num_gp,2) = j * dy
    838                    location(num_gp,3) = k * dz + 0.5 * dz
     875                   location(num_gp,3) = k * dz + 0.5_wp * dz
    839876                   ei(num_gp)     = e(k+1,j,i)
    840877                   dissi(num_gp)  = diss(k+1,j,i)
    841                    de_dxi(num_gp) = 0.0
     878                   de_dxi(num_gp) = 0.0_wp
    842879                   de_dyi(num_gp) = de_dy(k+1,j,i)
    843880                   de_dzi(num_gp) = de_dz(k+1,j,i)
     
    847884                     gp_outside_of_building(3) == 0 )  THEN
    848885                   num_gp = num_gp + 1
    849                    location(num_gp,1) = i * dx + 0.5 * dx
     886                   location(num_gp,1) = i * dx + 0.5_wp * dx
    850887                   location(num_gp,2) = j * dy
    851                    location(num_gp,3) = k * dz + 0.5 * dz
     888                   location(num_gp,3) = k * dz + 0.5_wp * dz
    852889                   ei(num_gp)     = e(k+1,j,i+1)
    853890                   dissi(num_gp)  = diss(k+1,j,i+1)
    854                    de_dxi(num_gp) = 0.0
     891                   de_dxi(num_gp) = 0.0_wp
    855892                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
    856893                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
     
    864901                   num_gp = num_gp + 1
    865902                   location(num_gp,1) = (i+1) * dx
    866                    location(num_gp,2) = j * dy + 0.5 * dy
    867                    location(num_gp,3) = k * dz + 0.5 * dz
     903                   location(num_gp,2) = j * dy + 0.5_wp * dy
     904                   location(num_gp,3) = k * dz + 0.5_wp * dz
    868905                   ei(num_gp)     = e(k+1,j,i+1)
    869906                   dissi(num_gp)  = diss(k+1,j,i+1)
    870907                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
    871                    de_dyi(num_gp) = 0.0
     908                   de_dyi(num_gp) = 0.0_wp
    872909                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
    873910                ENDIF
     
    877914                   num_gp = num_gp + 1
    878915                   location(num_gp,1) = (i+1) * dx
    879                    location(num_gp,2) = j * dy + 0.5 * dy
    880                    location(num_gp,3) = k * dz + 0.5 * dz
     916                   location(num_gp,2) = j * dy + 0.5_wp * dy
     917                   location(num_gp,3) = k * dz + 0.5_wp * dz
    881918                   ei(num_gp)     = e(k+1,j+1,i+1)
    882919                   dissi(num_gp)  = diss(k+1,j+1,i+1)
    883920                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    884                    de_dyi(num_gp) = 0.0
     921                   de_dyi(num_gp) = 0.0_wp
    885922                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    886923                ENDIF
     
    892929                     gp_outside_of_building(8) == 0 )  THEN
    893930                   num_gp = num_gp + 1
    894                    location(num_gp,1) = i * dx + 0.5 * dx
     931                   location(num_gp,1) = i * dx + 0.5_wp * dx
    895932                   location(num_gp,2) = (j+1) * dy
    896                    location(num_gp,3) = k * dz + 0.5 * dz
     933                   location(num_gp,3) = k * dz + 0.5_wp * dz
    897934                   ei(num_gp)     = e(k+1,j+1,i)
    898935                   dissi(num_gp)  = diss(k+1,j+1,i)
    899                    de_dxi(num_gp) = 0.0
     936                   de_dxi(num_gp) = 0.0_wp
    900937                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
    901938                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
     
    905942                     gp_outside_of_building(4) == 0 )  THEN
    906943                   num_gp = num_gp + 1
    907                    location(num_gp,1) = i * dx + 0.5 * dx
     944                   location(num_gp,1) = i * dx + 0.5_wp * dx
    908945                   location(num_gp,2) = (j+1) * dy
    909                    location(num_gp,3) = k * dz + 0.5 * dz
     946                   location(num_gp,3) = k * dz + 0.5_wp * dz
    910947                   ei(num_gp)     = e(k+1,j+1,i+1)
    911948                   dissi(num_gp)  = diss(k+1,j+1,i+1)
    912                    de_dxi(num_gp) = 0.0
     949                   de_dxi(num_gp) = 0.0_wp
    913950                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    914951                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     
    922959                   num_gp = num_gp + 1
    923960                   location(num_gp,1) = i * dx
    924                    location(num_gp,2) = j * dy + 0.5 * dy
    925                    location(num_gp,3) = k * dz + 0.5 * dz
     961                   location(num_gp,2) = j * dy + 0.5_wp * dy
     962                   location(num_gp,3) = k * dz + 0.5_wp * dz
    926963                   ei(num_gp)     = e(k+1,j,i)
    927964                   dissi(num_gp)  = diss(k+1,j,i)
    928965                   de_dxi(num_gp) = de_dx(k+1,j,i)
    929                    de_dyi(num_gp) = 0.0
     966                   de_dyi(num_gp) = 0.0_wp
    930967                   de_dzi(num_gp) = de_dz(k+1,j,i)
    931968                ENDIF
     
    935972                   num_gp = num_gp + 1
    936973                   location(num_gp,1) = i * dx
    937                    location(num_gp,2) = j * dy + 0.5 * dy
    938                    location(num_gp,3) = k * dz + 0.5 * dz
     974                   location(num_gp,2) = j * dy + 0.5_wp * dy
     975                   location(num_gp,3) = k * dz + 0.5_wp * dz
    939976                   ei(num_gp)     = e(k+1,j+1,i)
    940977                   dissi(num_gp)  = diss(k+1,j+1,i)
    941978                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
    942                    de_dyi(num_gp) = 0.0
     979                   de_dyi(num_gp) = 0.0_wp
    943980                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
    944981                ENDIF
     
    958995                   de_dxi(num_gp) = de_dx(k+1,j,i)
    959996                   de_dyi(num_gp) = de_dy(k+1,j,i)
    960                    de_dzi(num_gp) = 0.0
     997                   de_dzi(num_gp) = 0.0_wp
    961998                ENDIF
    962999
     
    9751012                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
    9761013                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
    977                    de_dzi(num_gp) = 0.0
     1014                   de_dzi(num_gp) = 0.0_wp
    9781015                ENDIF
    9791016
     
    9921029                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
    9931030                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
    994                    de_dzi(num_gp) = 0.0
     1031                   de_dzi(num_gp) = 0.0_wp
    9951032                ENDIF
    9961033
     
    10091046                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    10101047                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    1011                    de_dzi(num_gp) = 0.0
     1048                   de_dzi(num_gp) = 0.0_wp
    10121049                ENDIF
    10131050
     
    10191056!--                building, it follows that the values at the particle
    10201057!--                location are the same as the gridpoint values
    1021                    e_int     = ei(num_gp)
    1022                    diss_int  = dissi(num_gp)
    1023                    de_dx_int = de_dxi(num_gp)
    1024                    de_dy_int = de_dyi(num_gp)
    1025                    de_dz_int = de_dzi(num_gp)
     1058                   e_int(n)     = ei(num_gp)
     1059                   diss_int(n)  = dissi(num_gp)
     1060                   de_dx_int(n) = de_dxi(num_gp)
     1061                   de_dy_int(n) = de_dyi(num_gp)
     1062                   de_dz_int(n) = de_dzi(num_gp)
    10261063                ELSE IF ( num_gp > 1 )  THEN
    10271064
    1028                    d_sum = 0.0
     1065                   d_sum = 0.0_wp
    10291066!
    10301067!--                Evaluation of the distances between the gridpoints
     
    10341071                      d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
    10351072                                   + ( particles(n)%y-location(agp,2) )**2  &
    1036                                    + ( particles(n)%z-location(agp,3) )**2
     1073                                   + ( zv(n)-location(agp,3) )**2
    10371074                      d_sum        = d_sum + d_gp_pl(agp)
    10381075                   ENDDO
     
    10401077!
    10411078!--                Finally the interpolation can be carried out
    1042                    e_int     = 0.0
    1043                    diss_int  = 0.0
    1044                    de_dx_int = 0.0
    1045                    de_dy_int = 0.0
    1046                    de_dz_int = 0.0
     1079                   e_int(n)     = 0.0_wp
     1080                   diss_int(n)  = 0.0_wp
     1081                   de_dx_int(n) = 0.0_wp
     1082                   de_dy_int(n) = 0.0_wp
     1083                   de_dz_int(n) = 0.0_wp
    10471084                   DO  agp = 1, num_gp
    1048                       e_int     = e_int     + ( d_sum - d_gp_pl(agp) ) * &
     1085                      e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
    10491086                                             ei(agp) / ( (num_gp-1) * d_sum )
    1050                       diss_int  = diss_int  + ( d_sum - d_gp_pl(agp) ) * &
     1087                      diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
    10511088                                          dissi(agp) / ( (num_gp-1) * d_sum )
    1052                       de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
     1089                      de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    10531090                                         de_dxi(agp) / ( (num_gp-1) * d_sum )
    1054                       de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
     1091                      de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    10551092                                         de_dyi(agp) / ( (num_gp-1) * d_sum )
    1056                       de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
     1093                      de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    10571094                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
    10581095                   ENDDO
     
    10611098
    10621099             ENDIF
    1063 
    1064           ENDIF 
    1065 
    1066 !
    1067 !--       Vertically interpolate the horizontally averaged SGS TKE and
    1068 !--       resolved-scale velocity variances and use the interpolated values
    1069 !--       to calculate the coefficient fs, which is a measure of the ratio
    1070 !--       of the subgrid-scale turbulent kinetic energy to the total amount
    1071 !--       of turbulent kinetic energy.
    1072           IF ( k == 0 )  THEN
    1073              e_mean_int = hom(0,1,8,0)
    1074           ELSE
    1075              e_mean_int = hom(k,1,8,0) +                                    &
    1076                                         ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
    1077                                         ( zu(k+1) - zu(k) ) *               &
    1078                                         ( particles(n)%z - zu(k) )
    1079           ENDIF
    1080 
    1081           kw = particles(n)%z / dz
    1082 
    1083           IF ( k == 0 )  THEN
    1084              aa  = hom(k+1,1,30,0)  * ( particles(n)%z / &
    1085                                       ( 0.5 * ( zu(k+1) - zu(k) ) ) )
    1086              bb  = hom(k+1,1,31,0)  * ( particles(n)%z / &
    1087                                       ( 0.5 * ( zu(k+1) - zu(k) ) ) )
    1088              cc  = hom(kw+1,1,32,0) * ( particles(n)%z / &
    1089                                       ( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
    1090           ELSE
    1091              aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
    1092                         ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
    1093              bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
    1094                         ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
    1095              cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
    1096                         ( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
    1097           ENDIF
    1098 
    1099           vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
    1100 
    1101           fs_int = ( 2.0_wp / 3.0_wp ) * e_mean_int / &
    1102                       ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
    1103 
     1100          ENDDO
     1101       ENDIF
     1102
     1103       DO nb = 0,7
     1104          i = ip + block_offset(nb)%i_off
     1105          j = jp + block_offset(nb)%j_off
     1106          k = kp + block_offset(nb)%k_off
     1107
     1108          DO  n = start_index(nb), end_index(nb)
     1109!
     1110!--          Vertical interpolation of the horizontally averaged SGS TKE and
     1111!--          resolved-scale velocity variances and use the interpolated values
     1112!--          to calculate the coefficient fs, which is a measure of the ratio
     1113!--          of the subgrid-scale turbulent kinetic energy to the total amount
     1114!--          of turbulent kinetic energy.
     1115             IF ( k == 0 )  THEN
     1116                e_mean_int = hom(0,1,8,0)
     1117             ELSE
     1118                e_mean_int = hom(k,1,8,0) +                                    &
     1119                                           ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
     1120                                           ( zu(k+1) - zu(k) ) *               &
     1121                                           ( zv(n) - zu(k) )
     1122             ENDIF
     1123
     1124 !           kw = particles(n)%z / dz
     1125             kw = kp-1                                    ! ok for ocean??? ( + offset_ocean_nzt_m1 ???)
     1126
     1127             IF ( k == 0 )  THEN
     1128                aa  = hom(k+1,1,30,0)  * ( zv(n) / &
     1129                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
     1130                bb  = hom(k+1,1,31,0)  * ( zv(n) / &
     1131                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
     1132                cc  = hom(kw+1,1,32,0) * ( zv(n) / &
     1133                                         ( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
     1134             ELSE
     1135                aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) *    &
     1136                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
     1137                bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) *    &
     1138                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
     1139                cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *   &
     1140                           ( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
     1141             ENDIF
     1142
     1143             vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
     1144!
     1145!--          Needed to avoid NaN particle velocities. The value of 1.0 is just
     1146!--          an educated guess for the given case.
     1147             IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp )  THEN
     1148                fs_int(n) = 1.0_wp
     1149             ELSE
     1150                fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int /                 &
     1151                            ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
     1152             ENDIF
     1153
     1154          ENDDO
     1155       ENDDO
     1156
     1157       DO  n = 1, number_of_particles
     1158
     1159         rg(n,1) = random_gauss( iran_part, 5.0_wp )
     1160         rg(n,2) = random_gauss( iran_part, 5.0_wp )
     1161         rg(n,3) = random_gauss( iran_part, 5.0_wp )
     1162
     1163       ENDDO
     1164
     1165       DO  n = 1, number_of_particles
    11041166!
    11051167!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
    1106           lagr_timescale = ( 4.0 * e_int ) / &
    1107                            ( 3.0 * fs_int * c_0 * diss_int )
     1168          lagr_timescale = ( 4.0_wp * e_int(n) ) / &
     1169                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )
    11081170
    11091171!
     
    11111173!--       complete the current LES timestep.
    11121174          dt_gap = dt_3d - particles(n)%dt_sum
    1113           dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
     1175          dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
    11141176
    11151177!
    11161178!--       The particle timestep should not be too small in order to prevent
    11171179!--       the number of particle timesteps of getting too large
    1118           IF ( dt_particle < dt_min_part   .AND.  dt_min_part < dt_gap ) &
    1119           THEN
    1120              dt_particle = dt_min_part
     1180          IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
     1181             dt_particle(n) = dt_min_part
    11211182          ENDIF
    11221183
    11231184!
    11241185!--       Calculate the SGS velocity components
    1125           IF ( particles(n)%age == 0.0 )  THEN
     1186          IF ( particles(n)%age == 0.0_wp )  THEN
    11261187!
    11271188!--          For new particles the SGS components are derived from the SGS
     
    11291190!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
    11301191!--          from becoming unrealistically large.
    1131              particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * &
    1132                                         ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp )
    1133              particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * &
    1134                                         ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp )
    1135              particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * &
    1136                                         ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp )
     1192             particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wfu_part * e_int(n) ) *  &
     1193                                  ( rg(n,1) - 1.0_wp )
     1194             particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wfv_part * e_int(n) ) *  &
     1195                                  ( rg(n,2) - 1.0_wp )
     1196             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wfw_part * e_int(n) ) *  &
     1197                                  ( rg(n,3) - 1.0_wp )
    11371198
    11381199          ELSE
    1139 
    11401200!
    11411201!--          Restriction of the size of the new timestep: compared to the
     
    11431203
    11441204             dt_particle_m = particles(n)%age - particles(n)%age_m
    1145              IF ( dt_particle > 2.0 * dt_particle_m )  THEN
    1146                 dt_particle = 2.0 * dt_particle_m
     1205             IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
     1206                dt_particle(n) = 2.0_wp * dt_particle_m
    11471207             ENDIF
    11481208
     
    11531213!--          As negative values for the subgrid TKE are not allowed, the
    11541214!--          change of the subgrid TKE with time cannot be smaller than
    1155 !--          -e_int/dt_particle. This value is used as a lower boundary
     1215!--          -e_int(n)/dt_particle. This value is used as a lower boundary
    11561216!--          value for the change of TKE
    11571217
    1158              de_dt_min = - e_int / dt_particle
    1159 
    1160              de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
     1218             de_dt_min = - e_int(n) / dt_particle(n)
     1219
     1220             de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
    11611221
    11621222             IF ( de_dt < de_dt_min )  THEN
     
    11641224             ENDIF
    11651225
    1166              particles(n)%rvar1 = particles(n)%rvar1 - fs_int * c_0 *          &
    1167                                   diss_int * particles(n)%rvar1 * dt_particle /&
    1168                                   ( 4.0 * sgs_wfu_part * e_int ) +             &
    1169                                   ( 2.0 * sgs_wfu_part * de_dt *               &
    1170                                     particles(n)%rvar1 /                       &
    1171                                     ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
    1172                                   ) * dt_particle / 2.0          +             &
    1173                                   SQRT( fs_int * c_0 * diss_int ) *            &
    1174                                   ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp ) *   &
    1175                                   SQRT( dt_particle )
    1176 
    1177              particles(n)%rvar2 = particles(n)%rvar2 - fs_int * c_0 *          &
    1178                                   diss_int * particles(n)%rvar2 * dt_particle /&
    1179                                   ( 4.0 * sgs_wfv_part * e_int ) +             &
    1180                                   ( 2.0 * sgs_wfv_part * de_dt *               &
    1181                                     particles(n)%rvar2 /                       &
    1182                                     ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
    1183                                   ) * dt_particle / 2.0_wp          +             &
    1184                                   SQRT( fs_int * c_0 * diss_int ) *            &
    1185                                   ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp ) *   &
    1186                                   SQRT( dt_particle )
    1187 
    1188              particles(n)%rvar3 = particles(n)%rvar3 - fs_int * c_0 *          &
    1189                                   diss_int * particles(n)%rvar3 * dt_particle /&
    1190                                   ( 4.0 * sgs_wfw_part * e_int ) +             &
    1191                                   ( 2.0 * sgs_wfw_part * de_dt *               &
    1192                                     particles(n)%rvar3 /                       &
    1193                                     ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
    1194                                   ) * dt_particle / 2.0_wp         +             &
    1195                                   SQRT( fs_int * c_0 * diss_int ) *            &
    1196                                   ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp ) *   &
    1197                                   SQRT( dt_particle )
     1226             particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 *       &
     1227                           diss_int(n) * particles(n)%rvar1 * dt_particle(n) / &
     1228                           ( 4.0_wp * sgs_wfu_part * e_int(n) ) +              &
     1229                           ( 2.0_wp * sgs_wfu_part * de_dt *                   &
     1230                             particles(n)%rvar1 /                              &
     1231                             ( 2.0_wp * sgs_wfu_part * e_int(n) ) +            &
     1232                             de_dx_int(n)                                      &
     1233                           ) * dt_particle(n) / 2.0_wp          +              &
     1234                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
     1235                           ( rg(n,1) - 1.0_wp ) *                              &
     1236                           SQRT( dt_particle(n) )
     1237
     1238             particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 *       &
     1239                           diss_int(n) * particles(n)%rvar2 * dt_particle(n) / &
     1240                           ( 4.0_wp * sgs_wfv_part * e_int(n) ) +              &
     1241                           ( 2.0_wp * sgs_wfv_part * de_dt *                   &
     1242                             particles(n)%rvar2 /                              &
     1243                             ( 2.0_wp * sgs_wfv_part * e_int(n) ) +            &
     1244                             de_dy_int(n)                                      &
     1245                           ) * dt_particle(n) / 2.0_wp          +              &
     1246                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
     1247                           ( rg(n,2) - 1.0_wp ) *                              &
     1248                           SQRT( dt_particle(n) )
     1249
     1250             particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 *       &
     1251                           diss_int(n) * particles(n)%rvar3 * dt_particle(n) / &
     1252                           ( 4.0_wp * sgs_wfw_part * e_int(n) ) +              &
     1253                           ( 2.0_wp * sgs_wfw_part * de_dt *                   &
     1254                             particles(n)%rvar3 /                              &
     1255                             ( 2.0_wp * sgs_wfw_part * e_int(n) ) +            &
     1256                             de_dz_int(n)                                      &
     1257                           ) * dt_particle(n) / 2.0_wp          +              &
     1258                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
     1259                           ( rg(n,3) - 1.0_wp ) *                              &
     1260                           SQRT( dt_particle(n) )
    11981261
    11991262          ENDIF
    1200 
    1201           u_int = u_int + particles(n)%rvar1
    1202           v_int = v_int + particles(n)%rvar2
    1203           w_int = w_int + particles(n)%rvar3
     1263          u_int(n) = u_int(n) + particles(n)%rvar1
     1264          v_int(n) = v_int(n) + particles(n)%rvar2
     1265          w_int(n) = w_int(n) + particles(n)%rvar3
    12041266
    12051267!
    12061268!--       Store the SGS TKE of the current timelevel which is needed for
    12071269!--       for calculating the SGS particle velocities at the next timestep
    1208           particles(n)%e_m = e_int
    1209 
    1210        ELSE
    1211 !
    1212 !--       If no SGS velocities are used, only the particle timestep has to
    1213 !--       be set
    1214           dt_particle = dt_3d
    1215 
    1216        ENDIF
    1217 
    1218 !
    1219 !--    Store the old age of the particle ( needed to prevent that a
    1220 !--    particle crosses several PEs during one timestep, and for the
    1221 !--    evaluation of the subgrid particle velocity fluctuations )
    1222        particles(n)%age_m = particles(n)%age
    1223 
    1224 
    1225 !
    1226 !--    Particle advection
    1227        IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 )  THEN
    1228 !
    1229 !--       Pure passive transport (without particle inertia)
    1230           particles(n)%x = particles(n)%x + u_int * dt_particle
    1231           particles(n)%y = particles(n)%y + v_int * dt_particle
    1232           particles(n)%z = particles(n)%z + w_int * dt_particle
    1233 
    1234           particles(n)%speed_x = u_int
    1235           particles(n)%speed_y = v_int
    1236           particles(n)%speed_z = w_int
    1237 
    1238        ELSE
    1239 !
     1270          particles(n)%e_m = e_int(n)
     1271       ENDDO
     1272
     1273    ELSE
     1274!
     1275!--    If no SGS velocities are used, only the particle timestep has to
     1276!--    be set
     1277       dt_particle = dt_3d
     1278
     1279    ENDIF
     1280!
     1281!-- Store the old age of the particle ( needed to prevent that a
     1282!-- particle crosses several PEs during one timestep, and for the
     1283!-- evaluation of the subgrid particle velocity fluctuations )
     1284    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
     1285
     1286    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
     1287
     1288    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
     1289       DO  n = 1, number_of_particles
     1290
     1291!
     1292!--       Particle advection
     1293          IF ( dens_ratio(n) == 0.0_wp )  THEN
     1294!
     1295!--          Pure passive transport (without particle inertia)
     1296             particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
     1297             particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
     1298             particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
     1299
     1300             particles(n)%speed_x = u_int(n)
     1301             particles(n)%speed_y = v_int(n)
     1302             particles(n)%speed_z = w_int(n)
     1303
     1304          ELSE
     1305!
     1306!--          Transport of particles with inertia
     1307             particles(n)%x = particles(n)%x + particles(n)%speed_x * &
     1308                                               dt_particle(n)
     1309             particles(n)%y = particles(n)%y + particles(n)%speed_y * &
     1310                                               dt_particle(n)
     1311             particles(n)%z = particles(n)%z + particles(n)%speed_z * &
     1312                                               dt_particle(n)
     1313
     1314!
     1315!--          Update of the particle velocity
     1316             IF ( cloud_droplets )  THEN
     1317                exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /      &
     1318                           ( particles(n)%radius )**2 *                        &
     1319                           ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius &
     1320                             * SQRT( ( u_int(n) - particles(n)%speed_x )**2 +  &
     1321                                     ( v_int(n) - particles(n)%speed_y )**2 +  &
     1322                                     ( w_int(n) - particles(n)%speed_z )**2 )  &
     1323                                          / molecular_viscosity )**0.687_wp    &
     1324                           )
     1325
     1326                exp_term = EXP( -exp_arg * dt_particle(n) )
     1327             ELSEIF ( use_sgs_for_particles )  THEN
     1328                exp_arg  = particle_groups(particles(n)%group)%exp_arg
     1329                exp_term = EXP( -exp_arg * dt_particle(n) )
     1330             ELSE
     1331                exp_arg  = particle_groups(particles(n)%group)%exp_arg
     1332                exp_term = particle_groups(particles(n)%group)%exp_term
     1333             ENDIF
     1334             particles(n)%speed_x = particles(n)%speed_x * exp_term +          &
     1335                                    u_int(n) * ( 1.0_wp - exp_term )
     1336             particles(n)%speed_y = particles(n)%speed_y * exp_term +          &
     1337                                    v_int(n) * ( 1.0_wp - exp_term )
     1338             particles(n)%speed_z = particles(n)%speed_z * exp_term +          &
     1339                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) *  &
     1340                                    g / exp_arg ) * ( 1.0_wp - exp_term )
     1341          ENDIF
     1342
     1343       ENDDO
     1344   
     1345    ELSE
     1346
     1347       DO  n = 1, number_of_particles
     1348
    12401349!--       Transport of particles with inertia
    1241           particles(n)%x = particles(n)%x + particles(n)%speed_x * &
    1242                                             dt_particle
    1243           particles(n)%y = particles(n)%y + particles(n)%speed_y * &
    1244                                             dt_particle
    1245           particles(n)%z = particles(n)%z + particles(n)%speed_z * &
    1246                                             dt_particle
    1247 
     1350          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
     1351          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
     1352          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
    12481353!
    12491354!--       Update of the particle velocity
    1250           dens_ratio = particle_groups(particles(n)%group)%density_ratio
    12511355          IF ( cloud_droplets )  THEN
    1252              exp_arg  = 4.5 * dens_ratio * molecular_viscosity /        &
    1253                         ( particles(n)%radius )**2 *                    &
    1254                         ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius *    &
    1255                           SQRT( ( u_int - particles(n)%speed_x )**2 +   &
    1256                                 ( v_int - particles(n)%speed_y )**2 +   &
    1257                                 ( w_int - particles(n)%speed_z )**2 ) / &
    1258                                          molecular_viscosity )**0.687_wp   &
     1356
     1357             exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /         &
     1358                        ( particles(n)%radius )**2 *                           &
     1359                        ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius *  &
     1360                          SQRT( ( u_int(n) - particles(n)%speed_x )**2 +       &
     1361                                ( v_int(n) - particles(n)%speed_y )**2 +       &
     1362                                ( w_int(n) - particles(n)%speed_z )**2 ) /     &
     1363                                         molecular_viscosity )**0.687_wp       &
    12591364                        )
    1260              exp_term = EXP( -exp_arg * dt_particle )
     1365
     1366             exp_term = EXP( -exp_arg * dt_particle(n) )
    12611367          ELSEIF ( use_sgs_for_particles )  THEN
    12621368             exp_arg  = particle_groups(particles(n)%group)%exp_arg
    1263              exp_term = EXP( -exp_arg * dt_particle )
     1369             exp_term = EXP( -exp_arg * dt_particle(n) )
    12641370          ELSE
    12651371             exp_arg  = particle_groups(particles(n)%group)%exp_arg
     
    12671373          ENDIF
    12681374          particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
    1269                                  u_int * ( 1.0 - exp_term )
     1375                                 u_int(n) * ( 1.0_wp - exp_term )
    12701376          particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
    1271                                  v_int * ( 1.0 - exp_term )
     1377                                 v_int(n) * ( 1.0_wp - exp_term )
    12721378          particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
    1273                                  ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg )&
    1274                                  * ( 1.0 - exp_term )
    1275        ENDIF
    1276 
     1379                                 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
     1380                                 exp_arg ) * ( 1.0_wp - exp_term )
     1381       ENDDO
     1382
     1383    ENDIF
     1384
     1385    DO  n = 1, number_of_particles
    12771386!
    12781387!--    Increment the particle age and the total time that the particle
    12791388!--    has advanced within the particle timestep procedure
    1280        particles(n)%age    = particles(n)%age    + dt_particle
    1281        particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
     1389       particles(n)%age    = particles(n)%age    + dt_particle(n)
     1390       particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
    12821391
    12831392!
    12841393!--    Check whether there is still a particle that has not yet completed
    12851394!--    the total LES timestep
    1286        IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 )  THEN
     1395       IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
    12871396          dt_3d_reached_l = .FALSE.
    12881397       ENDIF
     
    12901399    ENDDO
    12911400
     1401    CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
    12921402
    12931403 END SUBROUTINE lpm_advec
  • palm/trunk/SOURCE/lpm_boundary_conds.f90

    r1321 r1359  
    2020! Current revisions:
    2121! -----------------
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2224!
    2325! Former revisions:
     
    6668
    6769    USE control_parameters,                                                    &
    68         ONLY:  dz, message_string, particle_maximum_age
     70        ONLY:  dz, message_string, particle_maximum_age, simulated_time
    6971
    7072    USE cpulog,                                                                &
     
    8183    USE particle_attributes,                                                   &
    8284        ONLY:  deleted_particles, deleted_tails, ibc_par_b, ibc_par_t,         &
    83                number_of_particles, particles, particle_mask,                  &
     85               number_of_particles, particles,                                 &
    8486               particle_tail_coordinates, particle_type, offset_ocean_nzt_m1,  &
    8587               tail_mask, use_particle_tails, use_sgs_for_particles
     
    158160
    159161          IF ( particles(n)%age > particle_maximum_age  .AND.  &
    160                particle_mask(n) )                              &
     162               particles(n)%particle_mask )                              &
    161163          THEN
    162              particle_mask(n)  = .FALSE.
     164             particles(n)%particle_mask  = .FALSE.
    163165             deleted_particles = deleted_particles + 1
    164166             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
     
    168170          ENDIF
    169171
    170           IF ( particles(n)%z >= zu(nz)  .AND.  particle_mask(n) )  THEN
     172          IF ( particles(n)%z >= zu(nz)  .AND.  particles(n)%particle_mask )  THEN
    171173             IF ( ibc_par_t == 1 )  THEN
    172174!
    173175!--             Particle absorption
    174                 particle_mask(n)  = .FALSE.
     176                particles(n)%particle_mask  = .FALSE.
    175177                deleted_particles = deleted_particles + 1
    176178                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
     
    181183!
    182184!--             Particle reflection
    183                 particles(n)%z       = 2.0 * zu(nz) - particles(n)%z
     185                particles(n)%z       = 2.0_wp * zu(nz) - particles(n)%z
    184186                particles(n)%speed_z = -particles(n)%speed_z
    185187                IF ( use_sgs_for_particles  .AND. &
    186                      particles(n)%rvar3 > 0.0 )  THEN
     188                     particles(n)%rvar3 > 0.0_wp )  THEN
    187189                   particles(n)%rvar3 = -particles(n)%rvar3
    188190                ENDIF
    189191                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    190                    particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
     192                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
    191193                                               particle_tail_coordinates(1,3,nn)
    192194                ENDIF
    193195             ENDIF
    194196          ENDIF
    195           IF ( particles(n)%z < zw(0)  .AND.  particle_mask(n) )  THEN
     197         
     198          IF ( particles(n)%z < zw(0)  .AND.  particles(n)%particle_mask )  THEN
    196199             IF ( ibc_par_b == 1 )  THEN
    197200!
    198201!--             Particle absorption
    199                 particle_mask(n)  = .FALSE.
     202                particles(n)%particle_mask  = .FALSE.
    200203                deleted_particles = deleted_particles + 1
    201204                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
     
    206209!
    207210!--             Particle reflection
    208                 particles(n)%z       = 2.0 * zw(0) - particles(n)%z
     211                particles(n)%z       = 2.0_wp * zw(0) - particles(n)%z
    209212                particles(n)%speed_z = -particles(n)%speed_z
    210213                IF ( use_sgs_for_particles  .AND. &
    211                      particles(n)%rvar3 < 0.0 )  THEN
     214                     particles(n)%rvar3 < 0.0_wp )  THEN
    212215                   particles(n)%rvar3 = -particles(n)%rvar3
    213216                ENDIF
    214217                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    215                    particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
     218                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
    216219                                               particle_tail_coordinates(1,3,nn)
    217220                ENDIF
    218221                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    219                    particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - &
     222                   particle_tail_coordinates(1,3,nn) = 2.0_wp * zw(0) - &
    220223                                               particle_tail_coordinates(1,3,nn)
    221224                ENDIF
     
    236239          dt_particle = particles(n)%age - particles(n)%age_m
    237240
    238           i2 = ( particles(n)%x + 0.5 * dx ) * ddx
    239           j2 = ( particles(n)%y + 0.5 * dy ) * ddy
     241          i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
     242          j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
    240243          k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
    241244
     
    251254             pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
    252255             pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
    253              i1 = ( pos_x_old + 0.5 * dx ) * ddx
    254              j1 = ( pos_y_old + 0.5 * dy ) * ddy
     256             i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
     257             j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
    255258             k1 = pos_z_old / dz + offset_ocean_nzt_m1
    256259
    257260!
    258261!--          Case 1
    259              IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )&
     262             IF ( particles(n)%x > pos_x_old  .AND. particles(n)%y > pos_y_old )&
    260263             THEN
    261264                t_index = 1
    262265
    263266                DO  i = i1, i2
    264                    xline      = i * dx + 0.5 * dx
     267                   xline      = i * dx + 0.5_wp * dx
    265268                   t(t_index) = ( xline - pos_x_old ) / &
    266269                                ( particles(n)%x - pos_x_old )
     
    269272
    270273                DO  j = j1, j2
    271                    yline      = j * dy + 0.5 * dy
     274                   yline      = j * dy + 0.5_wp * dy
    272275                   t(t_index) = ( yline - pos_y_old ) / &
    273276                                ( particles(n)%y - pos_y_old )
     
    314317                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    315318
    316                    i3 = ( pos_x + 0.5 * dx ) * ddx   
    317                    j3 = ( pos_y + 0.5 * dy ) * ddy
     319                   i3 = ( pos_x + 0.5_wp * dx ) * ddx   
     320                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
    318321                   k3 = pos_z / dz + offset_ocean_nzt_m1
    319322
     
    353356                      ENDIF
    354357
    355                       IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
     358                      IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
    356359                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    357360                         reflect_y = .TRUE.
     
    359362                      ENDIF
    360363
    361                       IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
     364                      IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
    362365                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    363366                         reflect_x = .TRUE.
     
    377380
    378381                DO  i = i1, i2
    379                    xline      = i * dx + 0.5 * dx
     382                   xline      = i * dx + 0.5_wp * dx
    380383                   t(t_index) = ( xline - pos_x_old ) / &
    381384                                ( particles(n)%x - pos_x_old )
     
    384387
    385388                DO  j = j1, j2, -1
    386                    yline      = j * dy - 0.5 * dy
     389                   yline      = j * dy - 0.5_wp * dy
    387390                   t(t_index) = ( pos_y_old - yline ) / &
    388391                                ( pos_y_old - particles(n)%y )
     
    428431                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    429432
    430                    i3 = ( pos_x + 0.5 * dx ) * ddx
    431                    j3 = ( pos_y + 0.5 * dy ) * ddy
     433                   i3 = ( pos_x + 0.5_wp * dx ) * ddx
     434                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
    432435                   k3 = pos_z / dz + offset_ocean_nzt_m1
    433436
     
    456459                      ENDIF
    457460
    458                       IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
     461                      IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
    459462                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    460463                         reflect_x = .TRUE.
     
    473476                      ENDIF
    474477
    475                       IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
     478                      IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
    476479                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
    477480                         reflect_y = .TRUE.
     
    491494
    492495                DO  i = i1, i2, -1
    493                    xline      = i * dx - 0.5 * dx
     496                   xline      = i * dx - 0.5_wp * dx
    494497                   t(t_index) = ( pos_x_old - xline ) / &
    495498                                ( pos_x_old - particles(n)%x )
     
    498501
    499502                DO  j = j1, j2
    500                    yline      = j * dy + 0.5 * dy
     503                   yline      = j * dy + 0.5_wp * dy
    501504                   t(t_index) = ( yline - pos_y_old ) / &
    502505                                ( particles(n)%y - pos_y_old )
     
    543546                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    544547
    545                    i3 = ( pos_x + 0.5 * dx ) * ddx
    546                    j3 = ( pos_y + 0.5 * dy ) * ddy
     548                   i3 = ( pos_x + 0.5_wp * dx ) * ddx
     549                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
    547550                   k3 = pos_z / dz + offset_ocean_nzt_m1
    548551
     
    571574                      ENDIF
    572575
    573                       IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
     576                      IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
    574577                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    575578                         reflect_y = .TRUE.
     
    588591                      ENDIF
    589592
    590                       IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
     593                      IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
    591594                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
    592595                         reflect_x = .TRUE.
     
    606609
    607610                DO  i = i1, i2, -1
    608                    xline      = i * dx - 0.5 * dx
     611                   xline      = i * dx - 0.5_wp * dx
    609612                   t(t_index) = ( pos_x_old - xline ) / &
    610613                                ( pos_x_old - particles(n)%x )
     
    613616
    614617                DO  j = j1, j2, -1
    615                    yline      = j * dy - 0.5 * dy
     618                   yline      = j * dy - 0.5_wp * dy
    616619                   t(t_index) = ( pos_y_old - yline ) / &
    617620                                ( pos_y_old - particles(n)%y )
     
    658661                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    659662
    660                    i3 = ( pos_x + 0.5 * dx ) * ddx   
    661                    j3 = ( pos_y + 0.5 * dy ) * ddy
     663                   i3 = ( pos_x + 0.5_wp * dx ) * ddx   
     664                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
    662665                   k3 = pos_z / dz + offset_ocean_nzt_m1
    663666
     
    686689                      ENDIF
    687690
    688                       IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
     691                      IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
    689692                           nzb_s_inner(j3,i5) /=0  .AND.          &
    690693                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
     
    704707                      ENDIF
    705708
    706                       IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
     709                      IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
    707710                           nzb_s_inner(j5,i3) /= 0  .AND.         &
    708711                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
     
    724727          IF ( reflect_z )  THEN
    725728
    726              particles(n)%z       = 2.0 * pos_z - prt_z
     729             particles(n)%z       = 2.0_wp * pos_z - prt_z
    727730             particles(n)%speed_z = - particles(n)%speed_z
    728731
     
    734737          ELSEIF ( reflect_y )  THEN
    735738
    736              particles(n)%y       = 2.0 * pos_y - prt_y
     739             particles(n)%y       = 2.0_wp * pos_y - prt_y
    737740             particles(n)%speed_y = - particles(n)%speed_y
    738741
     
    744747          ELSEIF ( reflect_x )  THEN
    745748
    746              particles(n)%x       = 2.0 * pos_x - prt_x
     749             particles(n)%x       = 2.0_wp * pos_x - prt_x
    747750             particles(n)%speed_x = - particles(n)%speed_x
    748751
  • palm/trunk/SOURCE/lpm_calc_liquid_water_content.f90

    r1321 r1359  
    2020! Current revisions:
    2121! ------------------
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2224!
    2325! Former revisions:
     
    6870
    6971    USE particle_attributes,                                                   &
    70         ONLY:  particles, prt_count, prt_start_index
     72        ONLY:  grid_particles, number_of_particles, particles, prt_count
    7173
    7274    IMPLICIT NONE
     
    7880    INTEGER(iwp) ::  psi !:
    7981
    80 
    8182    CALL cpu_log( log_point_s(45), 'lpm_calc_ql', 'start' )
    8283
    8384!
    8485!-- Set water content initially to zero
    85     ql = 0.0;  ql_v = 0.0;  ql_vp = 0.0
     86    ql = 0.0_wp;  ql_v = 0.0_wp;  ql_vp = 0.0_wp
    8687
    8788!
     
    8990    DO  i = nxl, nxr
    9091       DO  j = nys, nyn
    91           DO  k = nzb, nzt+1
     92          DO  k = nzb+1, nzt
     93
     94             number_of_particles = prt_count(k,j,i)
     95             IF ( number_of_particles <= 0 )  CYCLE
     96             particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    9297
    9398!
    9499!--          Calculate the total volume in the boxes (ql_v, weighting factor
    95100!--          has to beincluded)
    96              psi = prt_start_index(k,j,i)
    97              DO  n = psi, psi+prt_count(k,j,i)-1
     101             DO  n = 1, prt_count(k,j,i)
    98102                ql_v(k,j,i)  = ql_v(k,j,i)  + particles(n)%weight_factor *  &
    99103                                              particles(n)%radius**3
     
    102106!
    103107!--          Calculate the liquid water content
    104              IF ( ql_v(k,j,i) /= 0.0 )  THEN
    105                 ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi *           &
     108             IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
     109                ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333_wp * pi *           &
    106110                                        ql_v(k,j,i) /                       &
    107111                                        ( rho_surface * dx * dy * dz )
    108112
    109                 IF ( ql(k,j,i) < 0.0 ) THEN
     113                IF ( ql(k,j,i) < 0.0_wp ) THEN
    110114                   WRITE( message_string, * )  'LWC out of range: ' , &
    111                                                ql(k,j,i)
     115                                               ql(k,j,i),i,j,k
    112116                   CALL message( 'lpm_calc_liquid_water_content', '', 2, 2, &
    113117                                 -1, 6, 1 )
     
    116120             ELSE
    117121
    118                 ql(k,j,i) = 0.0
     122                ql(k,j,i) = 0.0_wp
    119123
    120124             ENDIF
  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r1347 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2324!
    2425! Former revisions:
     
    177178          rclass_lbound = LOG( 1.0E-6_wp )
    178179          rclass_ubound = LOG( 2.0E-4_wp )
    179           radclass(1)   = 1.0E-6
     180          radclass(1)   = 1.0E-6_wp
    180181          DO  i = 2, radius_classes
    181182             radclass(i) = EXP( rclass_lbound +                                &
    182                                 ( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /&
    183                                 ( radius_classes - 1.0 ) )
    184 !             IF ( myid == 0 )  THEN
    185 !                PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6
    186 !             ENDIF
     183                                ( rclass_ubound - rclass_lbound ) *            &
     184                                ( i - 1.0_wp ) / ( radius_classes - 1.0_wp ) )
    187185          ENDDO
    188186
     
    190188!--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
    191189          DO  i = 1, dissipation_classes
    192              epsclass(i) = 0.1 * REAL( i, KIND=wp ) / dissipation_classes
    193 !             IF ( myid == 0 )  THEN
    194 !                PRINT*, 'i=', i, ' eps = ', epsclass(i)
    195 !             ENDIF
     190             epsclass(i) = 0.1_wp * REAL( i, KIND=wp ) / dissipation_classes
    196191          ENDDO
    197192!
     
    205200
    206201             epsilon = epsclass(k)
    207              urms    = 2.02 * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
     202             urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
    208203
    209204             CALL turbsd
     
    240235
    241236             PRINT*, '*** Hall kernel'
    242              WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
     237             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6_wp, &
    243238                                              i = 1,radius_classes )
    244239             DO  j = 1, radius_classes
     
    250245                DO  i = 1, radius_classes
    251246                   DO  j = 1, radius_classes
    252                       IF ( hkernel(i,j) == 0.0 )  THEN
    253                          hwratio(i,j) = 9999999.9
     247                      IF ( hkernel(i,j) == 0.0_wp )  THEN
     248                         hwratio(i,j) = 9999999.9_wp
    254249                      ELSE
    255250                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
     
    259254
    260255                PRINT*, '*** epsilon = ', epsclass(k)
    261                 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
     256                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i) * 1.0E6_wp, &
    262257                                                 i = 1,radius_classes )
    263258                DO  j = 1, radius_classes
    264 !                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E6, &
    265 !                                       ( ckernel(i,j,k), i = 1,radius_classes )
    266                    WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j)*1.0E6, &
     259                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j) * 1.0E6_wp, &
    267260                                          ( hwratio(i,j), i = 1,radius_classes )
    268261                ENDDO
     
    292285
    293286       USE particle_attributes,                                                &
    294            ONLY:  prt_count, prt_start_index, radius_classes, wang_kernel
     287           ONLY:  prt_count, radius_classes, wang_kernel
    295288
    296289       IMPLICIT NONE
     
    305298
    306299
    307        pstart = prt_start_index(k1,j1,i1)
    308        pend   = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1
     300       pstart = 1
     301       pend   = prt_count(k1,j1,i1)
    309302       radius_classes = prt_count(k1,j1,i1)
    310303
     
    319312          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
    320313       ELSE
    321           epsilon = 0.0
     314          epsilon = 0.0_wp
    322315       ENDIF
    323        urms    = 2.02 * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
    324 
    325        IF ( wang_kernel  .AND.  epsilon > 1.0E-7 )  THEN
     316       urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
     317
     318       IF ( wang_kernel  .AND.  epsilon > 1.0E-7_wp )  THEN
    326319!
    327320!--       Call routines to calculate efficiencies for the Wang kernel
     
    442435       lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
    443436       tl        = urms**2 / epsilon                       ! in s
    444        lf        = 0.5 * urms**3 / epsilon                 ! in m
     437       lf        = 0.5_wp * urms**3 / epsilon              ! in m
    445438       tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
    446439       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp          ! in m
    447440       vk        = eta / tauk
    448441
    449        ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
    450        tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk   ! in s
     442       ao = ( 11.0_wp + 7.0_wp * lambda_re ) / ( 205.0_wp + lambda_re )
     443       tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk  ! in s
    451444
    452445       CALL fallg    ! gives winf in m/s
     
    461454       z   = tt / tl
    462455       be  = SQRT( 2.0_wp ) * lambda / lf
    463        bbb = SQRT( 1.0 - 2.0 * be**2 )
    464        d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
    465        e1  = lf * ( 1.0 + bbb ) * 0.5   ! in m
    466        d2  = ( 1.0 - bbb ) * 0.5 / bbb
    467        e2  = lf * ( 1.0 - bbb ) * 0.5   ! in m
    468        ccc = SQRT( 1.0 - 2.0 * z**2 )
    469        b1  = ( 1.0 + ccc ) * 0.5 / ccc
    470        c1  = tl * ( 1.0 + ccc ) * 0.5   ! in s
    471        b2  = ( 1.0 - ccc ) * 0.5 / ccc
    472        c2  = tl * ( 1.0 - ccc ) * 0.5   ! in s
     456       bbb = SQRT( 1.0_wp - 2.0_wp * be**2 )
     457       d1  = ( 1.0_wp + bbb ) / ( 2.0_wp * bbb )
     458       e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp   ! in m
     459       d2  = ( 1.0_wp - bbb ) * 0.5_wp / bbb
     460       e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp   ! in m
     461       ccc = SQRT( 1.0_wp - 2.0_wp * z**2 )
     462       b1  = ( 1.0_wp + ccc ) * 0.5_wp / ccc
     463       c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp   ! in s
     464       b2  = ( 1.0_wp - ccc ) * 0.5_wp / ccc
     465       c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp   ! in s
    473466
    474467       DO  i = 1, radius_classes
     
    509502                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
    510503             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
    511              v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in m**2/s**2
    512              wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy  ! in m**2/s**2
    513              IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
     504             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)   ! in m**2/s**2
     505             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy ! in m**2/s**2
     506             IF ( wrtur2xy < 0.0_wp )  wrtur2xy = 0.0_wp
    514507             wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
    515508             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) ) ! in m/s
     
    523516             ENDIF
    524517
    525              xx = -0.1988 * sst**4 + 1.5275 * sst**3 - 4.2942 * sst**2 + &
    526                    5.3406 * sst
    527              IF ( xx < 0.0 )  xx = 0.0
    528              yy = 0.1886 * EXP( 20.306_wp / lambda_re )
     518             xx = -0.1988_wp * sst**4 + 1.5275_wp * sst**3 - 4.2942_wp *      &
     519                   sst**2 + 5.3406_wp * sst
     520             IF ( xx < 0.0_wp )  xx = 0.0_wp
     521             yy = 0.1886_wp * EXP( 20.306_wp / lambda_re )
    529522
    530523             c1_gr  =  xx / ( g / vk * tauk )**yy
    531524
    532525             ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
    533              fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
     526             fao_gr = 20.115_wp * SQRT( ao_gr / lambda_re )
    534527             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
    535528
    536              grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 )
    537              IF ( grfin < 1.0 )  grfin = 1.0
    538 
    539              gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin           ! in cm**3/s
     529             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5_wp )
     530             IF ( grfin < 1.0_wp )  grfin = 1.0_wp
     531
     532             gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin        ! in cm**3/s
    540533             gck(j,i) = gck(i,j)
    541534
     
    559552       REAL(wp) ::  vsett !:
    560553
    561        aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
    562        phi_w = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
     554       aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b
     555       phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2  ! in s
    563556
    564557    END FUNCTION phi_w
     
    585578       REAL(wp) ::  vsett2 !:
    586579
    587        aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a
    588        aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a
    589        aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2
    590        aa4 = ( vsett2 / b )**2 - ( 1.0 / tau2 + 1.0 / a )**2
    591        aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a
    592        aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2
    593        zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 &
    594            + (4.0 / aa4 - 1.0 / aa5**2 - 1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6&
    595            + (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa2**2 + vsett2 / aa1**2 )&
    596            * 0.5 / b / aa3      ! in s**2
     580       aa1 = vsett2 / b - 1.0_wp / tau2 - 1.0_wp / a
     581       aa2 = vsett1 / b + 1.0_wp / tau1 + 1.0_wp / a
     582       aa3 = ( vsett1 - vsett2 ) / b + 1.0_wp / tau1 + 1.0_wp / tau2
     583       aa4 = ( vsett2 / b )**2 - ( 1.0_wp / tau2 + 1.0_wp / a )**2
     584       aa5 = vsett2 / b + 1.0_wp / tau2 + 1.0_wp / a
     585       aa6 = 1.0_wp / tau1 - 1.0_wp / a + ( 1.0_wp / tau2 + 1.0_wp / a) *      &
     586             vsett1 / vsett2
     587       zhi = (1.0_wp / aa1 - 1.0_wp / aa2 ) * ( vsett1 - vsett2 ) * 0.5_wp /   &
     588             b / aa3**2 + ( 4.0_wp / aa4 - 1.0_wp / aa5**2 - 1.0_wp / aa1**2 ) &
     589             * vsett2 * 0.5_wp / b /aa6 + ( 2.0_wp * ( b / aa2 - b / aa1 ) -   &
     590             vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3    ! in s**2
    597591
    598592    END FUNCTION zhi
     
    645639
    646640          first = .FALSE.
    647           b = (/  -0.318657E1,  0.992696E0, -0.153193E-2, -0.987059E-3, &
    648                  -0.578878E-3, 0.855176E-4, -0.327815E-5 /)
    649           c = (/  -0.500015E1,  0.523778E1,  -0.204914E1,   0.475294E0, &
    650                  -0.542819E-1, 0.238449E-2 /)
     641          b = (/  -0.318657E1_wp,   0.992696E0_wp,  -0.153193E-2_wp, &
     642                  -0.987059E-3_wp, -0.578878E-3_wp,  0.855176E-4_wp, &
     643                  -0.327815E-5_wp /)
     644          c = (/  -0.500015E1_wp,   0.523778E1_wp,  -0.204914E1_wp,   &
     645                   0.475294E0_wp,  -0.542819E-1_wp,  0.238449E-2_wp /)
    651646
    652647!
    653648!--       Parameter values for p = 1013,25 hPa and T = 293,15 K
    654           eta   = 1.818E-5         ! in kg/(m s)
    655           xlamb = 6.6E-8           ! in m
    656           rho_a = 1.204            ! in kg/m**3
    657           cunh  = 1.26 * xlamb     ! in m
    658           sigma = 0.07363          ! in kg/s**2
    659           stok  = 2.0  * g * ( rho_l - rho_a ) / ( 9.0 * eta ) ! in 1/(m s)
    660           stb   = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta)
     649          eta   = 1.818E-5_wp         ! in kg/(m s)
     650          xlamb = 6.6E-8_wp           ! in m
     651          rho_a = 1.204_wp            ! in kg/m**3
     652          cunh  = 1.26_wp * xlamb     ! in m
     653          sigma = 0.07363_wp          ! in kg/s**2
     654          stok  = 2.0_wp  * g * ( rho_l - rho_a ) / ( 9.0_wp * eta ) ! in 1/(m s)
     655          stb   = 32.0_wp * rho_a * ( rho_l - rho_a) * g / (3.0_wp * eta * eta)
    661656          phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
    662657          py    = phy**( 1.0_wp / 6.0_wp )
     
    666661       DO  j = 1, radius_classes
    667662
    668           IF ( radclass(j) <= 1.0E-5 ) THEN
     663          IF ( radclass(j) <= 1.0E-5_wp ) THEN
    669664
    670665             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) )
    671666
    672           ELSEIF ( radclass(j) > 1.0E-5  .AND.  radclass(j) <= 5.35E-4 )  THEN
     667          ELSEIF ( radclass(j) > 1.0E-5_wp  .AND.  radclass(j) <= 5.35E-4_wp )  THEN
    673668
    674669             x = LOG( stb * radclass(j)**3 )
    675              y = 0.0
     670             y = 0.0_wp
    676671
    677672             DO  i = 1, 7
     
    681676!--          Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418)
    682677!--          for correct version see (Beard, 1976)
    683              xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
    684 
    685              winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
    686 
    687           ELSEIF ( radclass(j) > 5.35E-4 )  THEN
    688 
    689              IF ( radclass(j) > 0.0035 )  THEN
    690                 bond = g * ( rho_l - rho_a ) * 0.0035**2 / sigma
     678             xrey = ( 1.0_wp + cunh / radclass(j) ) * EXP( y )
     679
     680             winf(j) = xrey * eta / ( 2.0_wp * rho_a * radclass(j) )
     681
     682          ELSEIF ( radclass(j) > 5.35E-4_wp )  THEN
     683
     684             IF ( radclass(j) > 0.0035_wp )  THEN
     685                bond = g * ( rho_l - rho_a ) * 0.0035_wp**2 / sigma
    691686             ELSE
    692687               bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma
    693688             ENDIF
    694689
    695              x = LOG( 16.0 * bond * py / 3.0_wp )
    696              y = 0.0
     690             x = LOG( 16.0_wp * bond * py / 3.0_wp )
     691             y = 0.0_wp
    697692
    698693             DO  i = 1, 6
     
    702697             xrey = py * EXP( y )
    703698
    704              IF ( radclass(j) > 0.0035 )  THEN
    705                 winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035_wp )
     699             IF ( radclass(j) > 0.0035_wp )  THEN
     700                winf(j) = xrey * eta / ( 2.0_wp * rho_a * 0.0035_wp )
    706701             ELSE
    707                 winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
     702                winf(j) = xrey * eta / ( 2.0_wp * rho_a * radclass(j) )
    708703             ENDIF
    709704
     
    752747
    753748         first = .FALSE.
    754          r0  = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., &
    755                   70.0, 100.0, 150.0, 200.0, 300.0 /)
    756          rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
    757                   0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, &
    758                   1.00 /)
    759 
    760          ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, &
    761                         0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/)
    762          ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, &
    763                         0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/)
    764          ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, &
    765                         0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/)
    766          ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, &
    767                         0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/)
    768          ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, &
    769                         0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/)
    770          ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, &
    771                         0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/)
    772          ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, &
    773                         0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/)
    774          ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, &
    775                         0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/)
    776          ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, &
    777                         0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/)
    778          ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, &
    779                         0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
    780          ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, &
    781                         0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
    782          ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, &
    783                         0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
    784          ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, &
    785                         0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
    786          ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, &
    787                         0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
    788          ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, &
    789                         0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/)
    790          ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, &
    791                         0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/)
    792          ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, &
    793                         0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/)
    794          ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, &
    795                         0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/)
    796          ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, &
    797                         0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/)
    798          ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, &
    799                         0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/)
    800          ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, &
    801                         1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/)
     749         r0  = (/   6.0_wp,   8.0_wp,  10.0_wp, 15.0_wp,  20.0_wp,  25.0_wp,  &
     750                   30.0_wp,  40.0_wp,  50.0_wp, 60.0_wp,  70.0_wp, 100.0_wp,  &
     751                  150.0_wp, 200.0_wp, 300.0_wp /)
     752
     753         rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp,       &
     754                  0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp,       &
     755                  0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp,       &
     756                  0.90_wp, 0.95_wp, 1.00_wp /)
     757
     758         ecoll(:,1)  = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, &
     759                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, &
     760                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp /)
     761         ecoll(:,2)  = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp, &
     762                          0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp, &
     763                          0.200_wp, 0.500_wp, 0.770_wp, 0.870_wp, 0.970_wp /)
     764         ecoll(:,3)  = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp, &
     765                          0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp, &
     766                          0.580_wp, 0.790_wp, 0.930_wp, 0.960_wp, 1.000_wp /)
     767         ecoll(:,4)  = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp, &
     768                          0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp, &
     769                          0.750_wp, 0.910_wp, 0.970_wp, 0.980_wp, 1.000_wp /)
     770         ecoll(:,5)  = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp, &
     771                          0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp, &
     772                          0.840_wp, 0.950_wp, 0.970_wp, 1.000_wp, 1.000_wp /)
     773         ecoll(:,6)  = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp, &
     774                          0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp, &
     775                          0.880_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     776         ecoll(:,7)  = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp, &
     777                          0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp, &
     778                          0.900_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     779         ecoll(:,8)  = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp, &
     780                          0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp, &
     781                          0.920_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     782         ecoll(:,9)  = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp, &
     783                          0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp, &
     784                          0.940_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     785         ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp, &
     786                          0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
     787                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     788         ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp, &
     789                          0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
     790                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     791         ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp, &
     792                          0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
     793                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     794         ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp, &
     795                          0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp, &
     796                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     797         ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp, &
     798                          0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp, &
     799                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     800         ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp, &
     801                          0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp, &
     802                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     803         ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp, &
     804                          0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp, &
     805                          0.970_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     806         ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp, &
     807                          0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp, &
     808                          1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
     809         ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp, &
     810                          0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp, &
     811                          1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp /)
     812         ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp, &
     813                          0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp, &
     814                          1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp /)
     815         ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, &
     816                          0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp, &
     817                          2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp /)
     818         ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, &
     819                          0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp, &
     820                          4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp /)
    802821       ENDIF
    803822
     
    832851                   pp = ( ( radclass(j) * 1.0E06_wp ) - r0(ir-1) ) / &
    833852                        ( r0(ir) - r0(ir-1) )
    834                    qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    835                    ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)  &
    836                              + pp * ( 1.0-qq ) * ecoll(ir,iq-1)          &
    837                              + qq * ( 1.0-pp ) * ecoll(ir-1,iq)          &
     853                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     854                   ec(j,i) = ( 1.0_wp - pp ) * ( 1.0_wp - qq )                 &
     855                             * ecoll(ir-1,iq-1)                                &
     856                             + pp * ( 1.0_wp - qq ) * ecoll(ir,iq-1)           &
     857                             + qq * ( 1.0_wp - pp ) * ecoll(ir-1,iq)           &
    838858                             + pp * qq * ecoll(ir,iq)
    839859                ELSE
    840860                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    841                    ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
     861                   ec(j,i) = ( 1.0_wp - qq ) * ecoll(1,iq-1) + qq * ecoll(1,iq)
    842862                ENDIF
    843863             ELSE
    844864                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    845                 ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
     865                ek = ( 1.0_wp - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
    846866                ec(j,i) = MIN( ek, 1.0_wp )
    847867             ENDIF
    848868
    849              IF ( ec(j,i) < 1.0E-20 )  ec(j,i) = 0.0
     869             IF ( ec(j,i) < 1.0E-20_wp )  ec(j,i) = 0.0_wp
    850870
    851871             ec(i,j) = ec(j,i)
     
    901921          first = .FALSE.
    902922
    903           r0  = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /)
    904           rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
     923          r0  = (/  10.0_wp, 20.0_wp, 30.0_wp, 40.0_wp, 50.0_wp, 60.0_wp,  &
     924                   100.0_wp /)
     925
     926          rat = (/ 0.0_wp, 0.1_wp, 0.2_wp, 0.3_wp, 0.4_wp, 0.5_wp, 0.6_wp, &
     927                   0.7_wp, 0.8_wp, 0.9_wp, 1.0_wp /)
    905928!
    906929!--       for 100 cm**2/s**3
    907           ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
    908           ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
    909           ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
    910           ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
    911           ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
    912           ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
    913           ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
    914           ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
    915           ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
    916           ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
    917           ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
     930          ecoll_100(:,1)  = (/  1.74_wp,   1.74_wp,   1.773_wp, 1.49_wp,  &
     931                                1.207_wp,  1.207_wp,  1.0_wp /)
     932          ecoll_100(:,2)  = (/  1.46_wp,   1.46_wp,   1.421_wp, 1.245_wp, &
     933                                1.069_wp,  1.069_wp,  1.0_wp /)
     934          ecoll_100(:,3)  = (/  1.32_wp,   1.32_wp,   1.245_wp, 1.123_wp, &
     935                                1.000_wp,  1.000_wp,  1.0_wp /)
     936          ecoll_100(:,4)  = (/  1.250_wp,  1.250_wp,  1.148_wp, 1.087_wp, &
     937                                1.025_wp,  1.025_wp,  1.0_wp /)
     938          ecoll_100(:,5)  = (/  1.186_wp,  1.186_wp,  1.066_wp, 1.060_wp, &
     939                                1.056_wp,  1.056_wp,  1.0_wp /)
     940          ecoll_100(:,6)  = (/  1.045_wp,  1.045_wp,  1.000_wp, 1.014_wp, &
     941                                1.028_wp,  1.028_wp,  1.0_wp /)
     942          ecoll_100(:,7)  = (/  1.070_wp,  1.070_wp,  1.030_wp, 1.038_wp, &
     943                                1.046_wp,  1.046_wp,  1.0_wp /)
     944          ecoll_100(:,8)  = (/  1.000_wp,  1.000_wp,  1.054_wp, 1.042_wp, &
     945                                1.029_wp,  1.029_wp,  1.0_wp /)
     946          ecoll_100(:,9)  = (/  1.223_wp,  1.223_wp,  1.117_wp, 1.069_wp, &
     947                                1.021_wp,  1.021_wp,  1.0_wp /)
     948          ecoll_100(:,10) = (/  1.570_wp,  1.570_wp,  1.244_wp, 1.166_wp, &
     949                                1.088_wp,  1.088_wp,  1.0_wp /)
     950          ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,  8.61_wp,  &
     951                                2.60_wp,   2.60_wp,   1.0_wp /)
    918952!
    919953!--       for 400 cm**2/s**3
    920           ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
    921           ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
    922           ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
    923           ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
    924           ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
    925           ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
    926           ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
    927           ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
    928           ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
    929           ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
    930           ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
     954          ecoll_400(:,1)  = (/  4.976_wp,  4.976_wp,  3.593_wp,  2.519_wp, &
     955                                1.445_wp,  1.445_wp,  1.0_wp /)
     956          ecoll_400(:,2)  = (/  2.984_wp,  2.984_wp,  2.181_wp,  1.691_wp, &
     957                                1.201_wp,  1.201_wp,  1.0_wp /)
     958          ecoll_400(:,3)  = (/  1.988_wp,  1.988_wp,  1.475_wp,  1.313_wp, &
     959                                1.150_wp,  1.150_wp,  1.0_wp /)
     960          ecoll_400(:,4)  = (/  1.490_wp,  1.490_wp,  1.187_wp,  1.156_wp, &
     961                                1.126_wp,  1.126_wp,  1.0_wp /)
     962          ecoll_400(:,5)  = (/  1.249_wp,  1.249_wp,  1.088_wp,  1.090_wp, &
     963                                1.092_wp,  1.092_wp,  1.0_wp /)
     964          ecoll_400(:,6)  = (/  1.139_wp,  1.139_wp,  1.130_wp,  1.091_wp, &
     965                                1.051_wp,  1.051_wp,  1.0_wp /)
     966          ecoll_400(:,7)  = (/  1.220_wp,  1.220_wp,  1.190_wp,  1.138_wp, &
     967                                1.086_wp,  1.086_wp,  1.0_wp /)
     968          ecoll_400(:,8)  = (/  1.325_wp,  1.325_wp,  1.267_wp,  1.165_wp, &
     969                                1.063_wp,  1.063_wp,  1.0_wp /)
     970          ecoll_400(:,9)  = (/  1.716_wp,  1.716_wp,  1.345_wp,  1.223_wp, &
     971                                1.100_wp,  1.100_wp,  1.0_wp /)
     972          ecoll_400(:,10) = (/  3.788_wp,  3.788_wp,  1.501_wp,  1.311_wp, &
     973                                1.120_wp,  1.120_wp,  1.0_wp /)
     974          ecoll_400(:,11) = (/ 36.52_wp,  36.52_wp,  19.16_wp,  22.80_wp,  &
     975                               26.0_wp,   26.0_wp,    1.0_wp /)
    931976
    932977       ENDIF
     
    9641009             ENDDO
    9651010
    966              y1 = 0.0001      ! for 0 m**2/s**3
     1011             y1 = 0.0001_wp      ! for 0 m**2/s**3
    9671012
    9681013             IF ( ir < 8 )  THEN
    9691014                IF ( ir >= 2 )  THEN
    970                    pp = ( radclass(j)*1.0E6 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
     1015                   pp = ( radclass(j)*1.0E6_wp - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
    9711016                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    972                    y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) + &
    973                                 pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   + &
    974                                 qq * ( 1.0-pp ) * ecoll_100(ir-1,iq)   + &
    975                                 pp * qq         * ecoll_100(ir,iq)
    976                    y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) + &
    977                                 pp * ( 1.0-qq ) * ecoll_400(ir,iq-1)   + &
    978                                 qq * ( 1.0-pp ) * ecoll_400(ir-1,iq)   + &
    979                                 pp * qq         * ecoll_400(ir,iq)
     1017                   y2 = ( 1.0_wp - pp ) * ( 1.0_wp - qq ) * ecoll_100(ir-1,iq-1) + &
     1018                                pp * ( 1.0_wp - qq ) * ecoll_100(ir,iq-1)        + &
     1019                                qq * ( 1.0_wp - pp ) * ecoll_100(ir-1,iq)        + &
     1020                                pp * qq              * ecoll_100(ir,iq)
     1021                   y3 = ( 1.0-pp ) * ( 1.0_wp - qq ) * ecoll_400(ir-1,iq-1)      + &
     1022                                pp * ( 1.0_wp - qq ) * ecoll_400(ir,iq-1)        + &
     1023                                qq * ( 1.0_wp - pp ) * ecoll_400(ir-1,iq)        + &
     1024                                pp * qq              * ecoll_400(ir,iq)
    9801025                ELSE
    9811026                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    982                    y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
    983                    y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
     1027                   y2 = ( 1.0_wp - qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
     1028                   y3 = ( 1.0_wp - qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
    9841029                ENDIF
    9851030             ELSE
    9861031                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    987                 y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
    988                 y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
     1032                y2 = ( 1.0_wp - qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
     1033                y3 = ( 1.0_wp - qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
    9891034             ENDIF
    9901035!
    9911036!--          Linear interpolation of dissipation rate in m**2/s**3
    992              IF ( epsilon <= 0.01 )  THEN
    993                 ecf(j,i) = ( epsilon - 0.01 ) / (   0.0 - 0.01 ) * y1 &
    994                          + ( epsilon -   0.0 ) / ( 0.01 -   0.0 ) * y2
    995              ELSEIF ( epsilon <= 0.06 )  THEN
    996                 ecf(j,i) = ( epsilon - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
    997                          + ( epsilon - 0.01 ) / ( 0.04 - 0.01 ) * y3
     1037             IF ( epsilon <= 0.01_wp )  THEN
     1038                ecf(j,i) = ( epsilon - 0.01_wp ) / ( 0.0_wp  - 0.01_wp ) * y1 &
     1039                         + ( epsilon - 0.0_wp  ) / ( 0.01_wp - 0.0_wp ) * y2
     1040             ELSEIF ( epsilon <= 0.06_wp )  THEN
     1041                ecf(j,i) = ( epsilon - 0.04_wp ) / ( 0.01_wp - 0.04_wp ) * y2 &
     1042                         + ( epsilon - 0.01_wp ) / ( 0.04_wp - 0.01_wp ) * y3
    9981043             ELSE
    999                 ecf(j,i) = (   0.06 - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
    1000                          + (   0.06 - 0.01 ) / ( 0.04 - 0.01 ) * y3
     1044                ecf(j,i) = ( 0.06_wp - 0.04_wp ) / ( 0.01_wp - 0.04_wp ) * y2 &
     1045                         + ( 0.06_wp - 0.01_wp ) / ( 0.04_wp - 0.01_wp ) * y3
    10011046             ENDIF
    10021047
    1003              IF ( ecf(j,i) < 1.0 )  ecf(j,i) = 1.0
     1048             IF ( ecf(j,i) < 1.0_wp )  ecf(j,i) = 1.0_wp
    10041049
    10051050             ecf(i,j) = ecf(j,i)
     
    10411086       REAL(wp)      ::  y       !:
    10421087 
    1043        REAL(wp), DIMENSION(1:9), SAVE      ::  collected_r = 0.0 !:
     1088       REAL(wp), DIMENSION(1:9), SAVE      ::  collected_r = 0.0_wp !:
    10441089       
    1045        REAL(wp), DIMENSION(1:19), SAVE     ::  collector_r = 0.0 !:
     1090       REAL(wp), DIMENSION(1:19), SAVE     ::  collector_r = 0.0_wp !:
    10461091       
    1047        REAL(wp), DIMENSION(1:9,1:19), SAVE ::  ef = 0.0          !:
     1092       REAL(wp), DIMENSION(1:9,1:19), SAVE ::  ef = 0.0_wp          !:
    10481093
    10491094       mean_rm = mean_r * 1.0E06_wp
     
    10521097       IF ( first )  THEN
    10531098
    1054           collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /)
    1055           collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0,  &
    1056                            150.0, 200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, &
    1057                            1400.0, 1800.0, 2400.0, 3000.0 /)
    1058 
    1059           ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0, &
    1060                       0.0 /)
    1061           ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12, 0.17, 0.17, 0.17, 0.0 /)
    1062           ef(:,3) = (/0.001, 0.001, 0.02,  0.13,  0.28, 0.37, 0.54, 0.55, 0.47/)
    1063           ef(:,4) = (/0.001, 0.001, 0.02,  0.23,  0.4,  0.55, 0.7,  0.75, 0.75/)
    1064           ef(:,5) = (/0.01,  0.01,  0.03,  0.3,   0.4,  0.58, 0.73, 0.75, 0.79/)
    1065           ef(:,6) = (/0.01,  0.01,  0.13,  0.38,  0.57, 0.68, 0.80, 0.86, 0.91/)
    1066           ef(:,7) = (/0.01,  0.085, 0.23,  0.52,  0.68, 0.76, 0.86, 0.92, 0.95/)
    1067           ef(:,8) = (/0.01,  0.14,  0.32,  0.60,  0.73, 0.81, 0.90, 0.94, 0.96/)
    1068           ef(:,9) = (/0.025, 0.25,  0.43,  0.66,  0.78, 0.83, 0.92, 0.95, 0.96/)
    1069           ef(:,10)= (/0.039, 0.3,   0.46,  0.69,  0.81, 0.87, 0.93, 0.95, 0.96/)
    1070           ef(:,11)= (/0.095, 0.33,  0.51,  0.72,  0.82, 0.87, 0.93, 0.96, 0.97/)
    1071           ef(:,12)= (/0.098, 0.36,  0.51,  0.73,  0.83, 0.88, 0.93, 0.96, 0.97/)
    1072           ef(:,13)= (/0.1,   0.36,  0.52,  0.74,  0.83, 0.88, 0.93, 0.96, 0.97/)
    1073           ef(:,14)= (/0.17,  0.4,   0.54,  0.72,  0.83, 0.88, 0.94, 0.98, 1.0 /)
    1074           ef(:,15)= (/0.15,  0.37,  0.52,  0.74,  0.82, 0.88, 0.94, 0.98, 1.0 /)
    1075           ef(:,16)= (/0.11,  0.34,  0.49,  0.71,  0.83, 0.88, 0.94, 0.95, 1.0 /)
    1076           ef(:,17)= (/0.08,  0.29,  0.45,  0.68,  0.8,  0.86, 0.96, 0.94, 1.0 /)
    1077           ef(:,18)= (/0.04,  0.22,  0.39,  0.62,  0.75, 0.83, 0.92, 0.96, 1.0 /)
    1078           ef(:,19)= (/0.02,  0.16,  0.33,  0.55,  0.71, 0.81, 0.90, 0.94, 1.0 /)
     1099          collected_r = (/    2.0_wp,    3.0_wp,    4.0_wp,    6.0_wp,    8.0_wp, &
     1100                             10.0_wp,   15.0_wp,   20.0_wp,   25.0_wp /)
     1101          collector_r = (/   10.0_wp,   20.0_wp,   30.0_wp,   40.0_wp,   50.0_wp, &
     1102                             60.0_wp,   80.0_wp,  100.0_wp,  150.0_wp,  200.0_wp, &
     1103                            300.0_wp,  400.0_wp,  500.0_wp,  600.0_wp, 1000.0_wp, &
     1104                           1400.0_wp, 1800.0_wp, 2400.0_wp, 3000.0_wp /)
     1105
     1106          ef(:,1)  = (/ 0.017_wp, 0.027_wp, 0.037_wp, 0.052_wp, 0.052_wp,      &
     1107                        0.052_wp, 0.052_wp, 0.0_wp,   0.0_wp /)
     1108          ef(:,2)  = (/ 0.001_wp, 0.016_wp, 0.027_wp, 0.060_wp, 0.12_wp,       &
     1109                        0.17_wp,  0.17_wp,  0.17_wp,  0.0_wp /)
     1110          ef(:,3)  = (/ 0.001_wp, 0.001_wp, 0.02_wp,  0.13_wp,  0.28_wp,       &
     1111                        0.37_wp,  0.54_wp,  0.55_wp,  0.47_wp/)
     1112          ef(:,4)  = (/ 0.001_wp, 0.001_wp, 0.02_wp,  0.23_wp,  0.4_wp,        &
     1113                        0.55_wp,  0.7_wp,   0.75_wp,  0.75_wp/)
     1114          ef(:,5)  = (/ 0.01_wp,  0.01_wp,  0.03_wp,  0.3_wp,   0.4_wp,        &
     1115                        0.58_wp,  0.73_wp,  0.75_wp,  0.79_wp/)
     1116          ef(:,6)  = (/ 0.01_wp,  0.01_wp,  0.13_wp,  0.38_wp,  0.57_wp,       &
     1117                        0.68_wp,  0.80_wp,  0.86_wp,  0.91_wp/)
     1118          ef(:,7)  = (/ 0.01_wp,  0.085_wp, 0.23_wp,  0.52_wp,  0.68_wp,       &
     1119                        0.76_wp,  0.86_wp,  0.92_wp,  0.95_wp/)
     1120          ef(:,8)  = (/ 0.01_wp,  0.14_wp,  0.32_wp,  0.60_wp,  0.73_wp,       &
     1121                        0.81_wp,  0.90_wp,  0.94_wp,  0.96_wp/)
     1122          ef(:,9)  = (/ 0.025_wp, 0.25_wp,  0.43_wp,  0.66_wp,  0.78_wp,       &
     1123                        0.83_wp,  0.92_wp,  0.95_wp,  0.96_wp/)
     1124          ef(:,10) = (/ 0.039_wp, 0.3_wp,   0.46_wp,  0.69_wp,  0.81_wp,       &
     1125                        0.87_wp,  0.93_wp,  0.95_wp,  0.96_wp/)
     1126          ef(:,11) = (/ 0.095_wp, 0.33_wp,  0.51_wp,  0.72_wp,  0.82_wp,       &
     1127                        0.87_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
     1128          ef(:,12) = (/ 0.098_wp, 0.36_wp,  0.51_wp,  0.73_wp,  0.83_wp,       &
     1129                        0.88_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
     1130          ef(:,13) = (/ 0.1_wp,   0.36_wp,  0.52_wp,  0.74_wp,  0.83_wp,       &
     1131                        0.88_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
     1132          ef(:,14) = (/ 0.17_wp,  0.4_wp,   0.54_wp,  0.72_wp,  0.83_wp,       &
     1133                        0.88_wp,  0.94_wp,  0.98_wp,  1.0_wp /)
     1134          ef(:,15) = (/ 0.15_wp,  0.37_wp,  0.52_wp,  0.74_wp,  0.82_wp,       &
     1135                        0.88_wp,  0.94_wp,  0.98_wp,  1.0_wp /)
     1136          ef(:,16) = (/ 0.11_wp,  0.34_wp,  0.49_wp,  0.71_wp,  0.83_wp,       &
     1137                        0.88_wp,  0.94_wp,  0.95_wp,  1.0_wp /)
     1138          ef(:,17) = (/ 0.08_wp,  0.29_wp,  0.45_wp,  0.68_wp,  0.8_wp,        &
     1139                        0.86_wp,  0.96_wp,  0.94_wp,  1.0_wp /)
     1140          ef(:,18) = (/ 0.04_wp,  0.22_wp,  0.39_wp,  0.62_wp,  0.75_wp,       &
     1141                        0.83_wp,  0.92_wp,  0.96_wp,  1.0_wp /)
     1142          ef(:,19) = (/ 0.02_wp,  0.16_wp,  0.33_wp,  0.55_wp,  0.71_wp,       &
     1143                        0.81_wp,  0.90_wp,  0.94_wp,  1.0_wp /)
    10791144
    10801145       ENDIF
     
    10881153       ENDDO
    10891154
    1090        IF ( rm < 10.0 )  THEN
    1091           e = 0.0
    1092        ELSEIF ( mean_rm < 2.0 )  THEN
    1093           e = 0.001
    1094        ELSEIF ( mean_rm >= 25.0 )  THEN
    1095           IF( j <= 2 )  e = 0.0
    1096           IF( j == 3 )  e = 0.47
    1097           IF( j == 4 )  e = 0.8
    1098           IF( j == 5 )  e = 0.9
    1099           IF( j >=6  )  e = 1.0
    1100        ELSEIF ( rm >= 3000.0 )  THEN
    1101           IF( i == 1 )  e = 0.02
    1102           IF( i == 2 )  e = 0.16
    1103           IF( i == 3 )  e = 0.33
    1104           IF( i == 4 )  e = 0.55
    1105           IF( i == 5 )  e = 0.71
    1106           IF( i == 6 )  e = 0.81
    1107           IF( i == 7 )  e = 0.90
    1108           IF( i >= 8 )  e = 0.94
     1155       IF ( rm < 10.0_wp )  THEN
     1156          e = 0.0_wp
     1157       ELSEIF ( mean_rm < 2.0_wp )  THEN
     1158          e = 0.001_wp
     1159       ELSEIF ( mean_rm >= 25.0_wp )  THEN
     1160          IF( j <= 2 )  e = 0.0_wp
     1161          IF( j == 3 )  e = 0.47_wp
     1162          IF( j == 4 )  e = 0.8_wp
     1163          IF( j == 5 )  e = 0.9_wp
     1164          IF( j >=6  )  e = 1.0_wp
     1165       ELSEIF ( rm >= 3000.0_wp )  THEN
     1166          IF( i == 1 )  e = 0.02_wp
     1167          IF( i == 2 )  e = 0.16_wp
     1168          IF( i == 3 )  e = 0.33_wp
     1169          IF( i == 4 )  e = 0.55_wp
     1170          IF( i == 5 )  e = 0.71_wp
     1171          IF( i == 6 )  e = 0.81_wp
     1172          IF( i == 7 )  e = 0.90_wp
     1173          IF( i >= 8 )  e = 0.94_wp
    11091174       ELSE
    11101175          x  = mean_rm - collected_r(i)
     
    11191184
    11201185          e = ( (gg-aa)*ef(i,j) + (gg-bb)*ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
    1121                 (gg-dd)*ef(i+1,j+1) ) / (3.0*gg)
     1186                (gg-dd)*ef(i+1,j+1) ) / (3.0_wp * gg)
    11221187       ENDIF
    11231188
  • palm/trunk/SOURCE/lpm_data_output_particles.f90

    r1329 r1359  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! netCDF output currently not available
     23! output of particle data in binary format adopted to new particle structure
    2324!
    2425! Former revisions:
     
    5253        ONLY:  cpu_log, log_point_s
    5354
     55    USE indices,                                                               &
     56        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     57
     58    USE kinds
     59
    5460    USE netcdf_control
    5561
    5662    USE particle_attributes,                                                   &
    57         ONLY:  maximum_number_of_particles, maximum_number_of_tailpoints,      &
    58                maximum_number_of_tails, number_of_particles, number_of_tails,  &
    59                particles, particle_tail_coordinates
     63        ONLY:  grid_particles, maximum_number_of_particles,                    &
     64               maximum_number_of_tailpoints, maximum_number_of_tails,          &
     65               number_of_particles, number_of_tails, particles,                &
     66               particle_tail_coordinates, prt_count
    6067
    6168    IMPLICIT NONE
    6269
     70    INTEGER(iwp) ::  ip !:
     71    INTEGER(iwp) ::  jp !:
     72    INTEGER(iwp) ::  kp !:
    6373
    6474    CALL cpu_log( log_point_s(40), 'lpm_data_output', 'start' )
     
    6979    CALL check_open( 85 )
    7080
    71     WRITE ( 85 )  simulated_time, maximum_number_of_particles, &
    72                   number_of_particles
    73     WRITE ( 85 )  particles
    74     WRITE ( 85 )  maximum_number_of_tailpoints, maximum_number_of_tails, &
    75                   number_of_tails
    76     IF ( maximum_number_of_tails > 0 )  THEN
    77        WRITE ( 85 )  particle_tail_coordinates, prt_time_count
    78     ENDIF
     81    WRITE ( 85 )  simulated_time
     82    WRITE ( 85 )  prt_count
     83         
     84    DO  ip = nxl, nxr
     85       DO  jp = nys, nyn
     86          DO  kp = nzb+1, nzt
     87             number_of_particles = prt_count(kp,jp,ip)
     88             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     89             IF ( number_of_particles <= 0 )  CYCLE
     90             WRITE ( 85 )  particles
     91          ENDDO
     92       ENDDO
     93    ENDDO
     94!
     95!-- particle tails currently not available
     96!    WRITE ( 85 )  maximum_number_of_tailpoints, maximum_number_of_tails, &
     97!                  number_of_tails
     98!    IF ( maximum_number_of_tails > 0 )  THEN
     99!       WRITE ( 85 )  particle_tail_coordinates, prt_time_count
     100!    ENDIF
    79101
    80102    CALL close_file( 85 )
     
    82104
    83105#if defined( __netcdf )
    84 !
    85 !-- Output in netCDF format
    86     CALL check_open( 108 )
    87 
    88 !
    89 !-- Update the NetCDF time axis
    90     prt_time_count = prt_time_count + 1
    91 
    92     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, &
    93                             (/ simulated_time /),        &
    94                             start = (/ prt_time_count /), count = (/ 1 /) )
    95     CALL handle_netcdf_error( 'lpm_data_output_particles', 1 )
    96 
    97 !
    98 !-- Output the real number of particles used
    99     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
    100                             (/ number_of_particles /),   &
    101                             start = (/ prt_time_count /), count = (/ 1 /) )
    102     CALL handle_netcdf_error( 'lpm_data_output_particles', 2 )
    103 
    104 !
    105 !-- Output all particle attributes
    106     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age,      &
    107                             start = (/ 1, prt_time_count /),               &
    108                             count = (/ maximum_number_of_particles /) )
    109     CALL handle_netcdf_error( 'lpm_data_output_particles', 3 )
    110 
    111     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize,&
    112                             start = (/ 1, prt_time_count /),                &
    113                             count = (/ maximum_number_of_particles /) )
    114     CALL handle_netcdf_error( 'lpm_data_output_particles', 4 )
    115 
    116     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, &
    117                             start = (/ 1, prt_time_count /),               &
    118                             count = (/ maximum_number_of_particles /) )
    119     CALL handle_netcdf_error( 'lpm_data_output_particles', 5 )
    120 
    121     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, &
    122                             start = (/ 1, prt_time_count /),               &
    123                             count = (/ maximum_number_of_particles /) )
    124     CALL handle_netcdf_error( 'lpm_data_output_particles', 6 )
    125 
    126     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, &
    127                             start = (/ 1, prt_time_count /),               &
    128                             count = (/ maximum_number_of_particles /) )
    129     CALL handle_netcdf_error( 'lpm_data_output_particles', 7 )
    130 
    131     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius,   &
    132                             start = (/ 1, prt_time_count /),               &
    133                             count = (/ maximum_number_of_particles /) )
    134     CALL handle_netcdf_error( 'lpm_data_output_particles', 8 )
    135 
    136     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x,  &
    137                             start = (/ 1, prt_time_count /),               &
    138                             count = (/ maximum_number_of_particles /) )
    139     CALL handle_netcdf_error( 'lpm_data_output_particles', 9 )
    140 
    141     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y,  &
    142                             start = (/ 1, prt_time_count /),               &
    143                             count = (/ maximum_number_of_particles /) )
    144     CALL handle_netcdf_error( 'lpm_data_output_particles', 10 )
    145 
    146     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z,  &
    147                             start = (/ 1, prt_time_count /),               &
    148                             count = (/ maximum_number_of_particles /) )
    149     CALL handle_netcdf_error( 'lpm_data_output_particles', 11 )
    150 
    151     nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),                     &
    152                             particles%weight_factor,                       &
    153                             start = (/ 1, prt_time_count /),               &
    154                             count = (/ maximum_number_of_particles /) )
    155     CALL handle_netcdf_error( 'lpm_data_output_particles', 12 )
    156 
    157     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x,       &
    158                             start = (/ 1, prt_time_count /),               &
    159                             count = (/ maximum_number_of_particles /) )
    160     CALL handle_netcdf_error( 'lpm_data_output_particles', 13 )
    161 
    162     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y,       &
    163                             start = (/ 1, prt_time_count /),               &
    164                             count = (/ maximum_number_of_particles /) )
    165     CALL handle_netcdf_error( 'lpm_data_output_particles', 14 )
    166 
    167     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z,       &
    168                             start = (/ 1, prt_time_count /),               &
    169                             count = (/ maximum_number_of_particles /) )
    170     CALL handle_netcdf_error( 'lpm_data_output_particles', 15 )
    171 
    172     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class,   &
    173                             start = (/ 1, prt_time_count /),               &
    174                             count = (/ maximum_number_of_particles /) )
    175     CALL handle_netcdf_error( 'lpm_data_output_particles', 16 )
    176 
    177     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group,   &
    178                             start = (/ 1, prt_time_count /),               &
    179                             count = (/ maximum_number_of_particles /) )
    180     CALL handle_netcdf_error( 'lpm_data_output_particles', 17 )
    181 
    182     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16),                    &
    183                             particles%tailpoints,                          &
    184                             start = (/ 1, prt_time_count /),               &
    185                             count = (/ maximum_number_of_particles /) )
    186     CALL handle_netcdf_error( 'lpm_data_output_particles', 18 )
    187 
    188     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, &
    189                             start = (/ 1, prt_time_count /),               &
    190                             count = (/ maximum_number_of_particles /) )
    191     CALL handle_netcdf_error( 'lpm_data_output_particles', 19 )
    192 
     106! !
     107! !-- Output in netCDF format
     108!     CALL check_open( 108 )
     109!
     110! !
     111! !-- Update the NetCDF time axis
     112!     prt_time_count = prt_time_count + 1
     113!
     114!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, &
     115!                             (/ simulated_time /),        &
     116!                             start = (/ prt_time_count /), count = (/ 1 /) )
     117!     CALL handle_netcdf_error( 'lpm_data_output_particles', 1 )
     118!
     119! !
     120! !-- Output the real number of particles used
     121!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
     122!                             (/ number_of_particles /),   &
     123!                             start = (/ prt_time_count /), count = (/ 1 /) )
     124!     CALL handle_netcdf_error( 'lpm_data_output_particles', 2 )
     125!
     126! !
     127! !-- Output all particle attributes
     128!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age,      &
     129!                             start = (/ 1, prt_time_count /),               &
     130!                             count = (/ maximum_number_of_particles /) )
     131!     CALL handle_netcdf_error( 'lpm_data_output_particles', 3 )
     132!
     133!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize,&
     134!                             start = (/ 1, prt_time_count /),                &
     135!                             count = (/ maximum_number_of_particles /) )
     136!     CALL handle_netcdf_error( 'lpm_data_output_particles', 4 )
     137!
     138!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, &
     139!                             start = (/ 1, prt_time_count /),               &
     140!                             count = (/ maximum_number_of_particles /) )
     141!     CALL handle_netcdf_error( 'lpm_data_output_particles', 5 )
     142!
     143!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, &
     144!                             start = (/ 1, prt_time_count /),               &
     145!                             count = (/ maximum_number_of_particles /) )
     146!     CALL handle_netcdf_error( 'lpm_data_output_particles', 6 )
     147!
     148!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, &
     149!                             start = (/ 1, prt_time_count /),               &
     150!                             count = (/ maximum_number_of_particles /) )
     151!     CALL handle_netcdf_error( 'lpm_data_output_particles', 7 )
     152!
     153!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius,   &
     154!                             start = (/ 1, prt_time_count /),               &
     155!                             count = (/ maximum_number_of_particles /) )
     156!     CALL handle_netcdf_error( 'lpm_data_output_particles', 8 )
     157!
     158!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x,  &
     159!                             start = (/ 1, prt_time_count /),               &
     160!                             count = (/ maximum_number_of_particles /) )
     161!     CALL handle_netcdf_error( 'lpm_data_output_particles', 9 )
     162!
     163!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y,  &
     164!                             start = (/ 1, prt_time_count /),               &
     165!                             count = (/ maximum_number_of_particles /) )
     166!     CALL handle_netcdf_error( 'lpm_data_output_particles', 10 )
     167!
     168!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z,  &
     169!                             start = (/ 1, prt_time_count /),               &
     170!                             count = (/ maximum_number_of_particles /) )
     171!     CALL handle_netcdf_error( 'lpm_data_output_particles', 11 )
     172!
     173!     nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),                     &
     174!                             particles%weight_factor,                       &
     175!                             start = (/ 1, prt_time_count /),               &
     176!                             count = (/ maximum_number_of_particles /) )
     177!     CALL handle_netcdf_error( 'lpm_data_output_particles', 12 )
     178!
     179!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x,       &
     180!                             start = (/ 1, prt_time_count /),               &
     181!                             count = (/ maximum_number_of_particles /) )
     182!     CALL handle_netcdf_error( 'lpm_data_output_particles', 13 )
     183!
     184!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y,       &
     185!                             start = (/ 1, prt_time_count /),               &
     186!                             count = (/ maximum_number_of_particles /) )
     187!     CALL handle_netcdf_error( 'lpm_data_output_particles', 14 )
     188!
     189!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z,       &
     190!                             start = (/ 1, prt_time_count /),               &
     191!                             count = (/ maximum_number_of_particles /) )
     192!     CALL handle_netcdf_error( 'lpm_data_output_particles', 15 )
     193!
     194!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class,   &
     195!                             start = (/ 1, prt_time_count /),               &
     196!                             count = (/ maximum_number_of_particles /) )
     197!     CALL handle_netcdf_error( 'lpm_data_output_particles', 16 )
     198!
     199!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group,   &
     200!                             start = (/ 1, prt_time_count /),               &
     201!                             count = (/ maximum_number_of_particles /) )
     202!     CALL handle_netcdf_error( 'lpm_data_output_particles', 17 )
     203!
     204!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16),                    &
     205!                             particles%tailpoints,                          &
     206!                             start = (/ 1, prt_time_count /),               &
     207!                             count = (/ maximum_number_of_particles /) )
     208!     CALL handle_netcdf_error( 'lpm_data_output_particles', 18 )
     209!
     210!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, &
     211!                             start = (/ 1, prt_time_count /),               &
     212!                             count = (/ maximum_number_of_particles /) )
     213!     CALL handle_netcdf_error( 'lpm_data_output_particles', 19 )
     214!
    193215#endif
    194216
  • palm/trunk/SOURCE/lpm_droplet_collision.f90

    r1323 r1359  
    1  SUBROUTINE lpm_droplet_collision
     1 SUBROUTINE lpm_droplet_collision (i,j,k)
    22
    33!--------------------------------------------------------------------------------!
     
    2020! Current revisions:
    2121! ------------------
    22 !
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2324!
    2425! Former revisions:
     
    7576!------------------------------------------------------------------------------!
    7677
     78
    7779    USE arrays_3d,                                                             &
    7880        ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
     
    103105    USE particle_attributes,                                                   &
    104106        ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
    105                palm_kernel, particles, particle_mask, particle_type,           &
    106                prt_count, prt_start_index, use_kernel_tables, wang_kernel
     107               palm_kernel, particles, particle_type,                          &
     108               prt_count, use_kernel_tables, wang_kernel
     109
     110    USE pegrid
    107111
    108112    IMPLICIT NONE
     
    124128    INTEGER(iwp) ::  rclass_s !:
    125129
     130    INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !:
     131
     132    LOGICAL, SAVE ::  first_flag = .TRUE. !:
     133
     134    TYPE(particle_type) :: tmp_particle   !:
     135
    126136    REAL(wp) ::  aa       !:
     137    REAL(wp) ::  auxn     !: temporary variables
     138    REAL(wp) ::  auxs     !: temporary variables
    127139    REAL(wp) ::  bb       !:
    128140    REAL(wp) ::  cc       !:
     
    158170    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !:
    159171
    160 
    161     TYPE(particle_type) ::  tmp_particle           !:
    162 
    163 
     172    REAL, DIMENSION(prt_count(k,j,i))    :: ck
     173    REAL, DIMENSION(prt_count(k,j,i))    :: r3v
     174    REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
     175    REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
    164176
    165177    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
    166178
    167     DO  i = nxl, nxr
    168        DO  j = nys, nyn
    169           DO  k = nzb+1, nzt
    170 !
    171 !--          Collision requires at least two particles in the box
    172              IF ( prt_count(k,j,i) > 1 )  THEN
    173 !
    174 !--             First, sort particles within the gridbox by their size,
    175 !--             using Shell's method (see Numerical Recipes)
    176 !--             NOTE: In case of using particle tails, the re-sorting of
    177 !--             ----  tails would have to be included here!
    178                 psi = prt_start_index(k,j,i) - 1
    179                 inc = 1
    180                 DO WHILE ( inc <= prt_count(k,j,i) )
    181                    inc = 3 * inc + 1
     179!
     180!-- Collision requires at least two particles in the box
     181    IF ( prt_count(k,j,i) > 1 )  THEN
     182!
     183!--    First, sort particles within the gridbox by their size,