Changeset 1822


Ignore:
Timestamp:
Apr 7, 2016 7:49:42 AM (5 years ago)
Author:
hoffmann
Message:

changes in LPM and bulk cloud microphysics

Location:
palm/trunk/SOURCE
Files:
5 deleted
40 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1818 r1822  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# Tails removed. lpm_release_set removed. calc_precipitation, impact_of_latent_heat
     23# removed.
    2324#
    2425# Former revisions:
     
    246247        advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \
    247248        buoyancy.f90 calc_liquid_water_content.f90 calc_mean_profile.f90 \
    248         calc_precipitation.f90 calc_radiation.f90 \
     249        calc_radiation.f90 \
    249250        check_for_restart.f90 check_open.f90 check_parameters.f90 \
    250251        close_file.f90 compute_vpt.f90 coriolis.f90 cpulog.f90 \
     
    257258        disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
    258259        exchange_horiz_2d.f90 fft_xy.f90 flow_statistics.f90 \
    259         global_min_max.f90 header.f90 impact_of_latent_heat.f90 \
     260        global_min_max.f90 header.f90 \
    260261        inflow_turbulence.f90 init_1d_model.f90 init_3d_model.f90 \
    261262        init_advec.f90 init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 \
     
    268269        lpm_data_output_particles.f90 lpm_droplet_collision.f90 \
    269270        lpm_droplet_condensation.f90 lpm_exchange_horiz.f90 \
    270         lpm_extend_tails.f90 \
    271         lpm_extend_tail_array.f90 lpm_init.f90 lpm_init_sgs_tke.f90 \
    272         lpm_pack_arrays.f90 lpm_read_restart_file.f90 lpm_release_set.f90 \
     271        lpm_init.f90 lpm_init_sgs_tke.f90 \
     272        lpm_pack_arrays.f90 lpm_read_restart_file.f90 \
    273273        lpm_set_attributes.f90 \
    274274        lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \
     
    349349calc_mean_profile.o: modules.o mod_kinds.o
    350350calc_liquid_water_content.o: modules.o mod_kinds.o
    351 calc_precipitation.o: modules.o mod_kinds.o
    352351calc_radiation.o: modules.o mod_kinds.o
    353352check_for_restart.o: modules.o mod_kinds.o pmc_interface.o
     
    395394   plant_canopy_model.o pmc_handle_communicator.o pmc_interface.o \
    396395   radiation_model.o spectrum.o subsidence.o
    397 impact_of_latent_heat.o: modules.o mod_kinds.o
    398396inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
    399397init_1d_model.o: modules.o mod_kinds.o
     
    418416local_tremain.o: modules.o cpulog.o mod_kinds.o
    419417local_tremain_ini.o: modules.o cpulog.o mod_kinds.o
    420 lpm.o: modules.o cpulog.o lpm_exchange_horiz.o lpm_pack_arrays.o mod_kinds.o \
    421         mod_particle_attributes.o
     418lpm.o: modules.o cpulog.o lpm_exchange_horiz.o lpm_init.o lpm_pack_arrays.o \
     419        mod_kinds.o mod_particle_attributes.o
    422420lpm_advec.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
    423421lpm_boundary_conds.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
     
    434432lpm_exchange_horiz.o: modules.o cpulog.o lpm_pack_arrays.o mod_kinds.o \
    435433        mod_particle_attributes.o netcdf_interface.o
    436 lpm_extend_tails.o: modules.o mod_kinds.o mod_particle_attributes.o
    437 lpm_extend_tail_array.o: modules.o mod_kinds.o mod_particle_attributes.o
    438434lpm_init.o: modules.o lpm_collision_kernels.o mod_kinds.o \
    439435        netcdf_interface.o random_function.o mod_particle_attributes.o \
     
    443439lpm_read_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o \
    444440        lpm_pack_arrays.o
    445 lpm_release_set.o: modules.o mod_kinds.o random_function.o \
    446         mod_particle_attributes.o lpm_init.o  random_generator_parallel.o
    447441lpm_set_attributes.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
    448442lpm_write_exchange_statistics.o: modules.o mod_kinds.o mod_particle_attributes.o
     
    480474prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_s_bc.o advec_u_pw.o \
    481475        advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o \
    482         advec_ws.o buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \
     476        advec_ws.o buoyancy.o calc_radiation.o coriolis.o \
    483477        cpulog.o diffusion_e.o diffusion_s.o diffusion_u.o diffusion_v.o diffusion_w.o \
    484         eqn_state_seawater.o impact_of_latent_heat.o mod_kinds.o microphysics.o \
     478        eqn_state_seawater.o mod_kinds.o microphysics.o \
    485479        nudging.o plant_canopy_model.o production_e.o radiation_model.o \
    486480        subsidence.o user_actions.o
  • palm/trunk/SOURCE/advec_ws.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! icloud_scheme removed, microphysics_seifert added
    2222!
    2323! Former revisions:
     
    222222
    223223       USE control_parameters,                                                 &
    224            ONLY:  cloud_physics, humidity, icloud_scheme, loop_optimization,   &
    225                   monotonic_adjustment, passive_scalar, precipitation, ocean,  &
    226                   ws_scheme_mom, ws_scheme_sca
     224           ONLY:  cloud_physics, humidity, loop_optimization,                  &
     225                  monotonic_adjustment, passive_scalar, microphysics_seifert,  &
     226                  ocean, ws_scheme_mom, ws_scheme_sca
    227227
    228228       USE indices,                                                            &
     
    274274          ENDIF
    275275
    276           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    277                precipitation )  THEN
     276          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    278277             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    279278             ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     
    331330             ENDIF
    332331
    333              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    334                   precipitation )  THEN
     332             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    335333                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
    336334                          diss_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
     
    365363   
    366364       USE control_parameters,                                                 &
    367            ONLY:  cloud_physics, humidity, icloud_scheme, passive_scalar,      &
    368                   precipitation, ocean, ws_scheme_mom, ws_scheme_sca
     365           ONLY:  cloud_physics, humidity, passive_scalar, ocean,              &
     366                  microphysics_seifert, ws_scheme_mom, ws_scheme_sca
    369367
    370368       USE kinds
     
    391389          sums_wspts_ws_l = 0.0_wp
    392390          IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0_wp
    393           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    394                precipitation )  THEN
     391          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    395392             sums_wsqrs_ws_l = 0.0_wp
    396393             sums_wsnrs_ws_l = 0.0_wp
  • palm/trunk/SOURCE/boundary_conds.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! icloud_scheme removed. microphyisics_seifert added.
    2222!
    2323! Former revisions:
     
    139139               cloud_physics, dt_3d, humidity,                                 &
    140140               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b,       &
    141                ibc_uv_t, icloud_scheme, inflow_l, inflow_n, inflow_r, inflow_s,&
    142                intermediate_timestep_count, large_scale_forcing, nest_domain,  &
    143                nest_bound_l, nest_bound_s, nudging, ocean,                     &
    144                outflow_l, outflow_n, outflow_r, outflow_s, passive_scalar,     &
    145                precipitation, tsc, use_cmax
     141               ibc_uv_t, inflow_l, inflow_n, inflow_r, inflow_s,               &
     142               intermediate_timestep_count, large_scale_forcing,               &
     143               microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
     144               nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
     145               passive_scalar, tsc, use_cmax
    146146
    147147    USE grid_variables,                                                        &
     
    319319       ENDIF
    320320
    321        IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
     321       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    322322!             
    323323!--       Surface conditions rain water (Dirichlet)
     
    373373       IF ( humidity  .OR.  passive_scalar )  THEN
    374374          q_p(:,nys-1,:) = q_p(:,nys,:)
    375           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    376                precipitation)  THEN
     375          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    377376             qr_p(:,nys-1,:) = qr_p(:,nys,:)
    378377             nr_p(:,nys-1,:) = nr_p(:,nys,:)
     
    384383       IF ( humidity  .OR.  passive_scalar )  THEN
    385384          q_p(:,nyn+1,:) = q_p(:,nyn,:)
    386           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    387                precipitation )  THEN
     385          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    388386             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
    389387             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
     
    395393       IF ( humidity  .OR.  passive_scalar )  THEN
    396394          q_p(:,:,nxl-1) = q_p(:,:,nxl)
    397           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    398                precipitation )  THEN
     395          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    399396             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
    400397             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
     
    406403       IF ( humidity .OR. passive_scalar )  THEN
    407404          q_p(:,:,nxr+1) = q_p(:,:,nxr)
    408           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
     405          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    409406             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
    410407             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
  • palm/trunk/SOURCE/calc_liquid_water_content.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! icloud_scheme removed. microphysics_seifert added.
    2222!
    2323! Former revisions:
     
    7373
    7474    USE control_parameters,                                                    &
    75         ONLY:  icloud_scheme, precipitation
     75        ONLY:  microphysics_seifert
    7676
    7777    USE indices,                                                               &
     
    122122!
    123123!--          Compute the liquid water content
    124              IF ( icloud_scheme == 0  .AND.  precipitation)  THEN
     124             IF ( microphysics_seifert )  THEN
    125125                IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN
    126126                   qc(k,j,i) = q(k,j,i) - q_s - qr(k,j,i)
     
    131131                   ql(k,j,i) = qr(k,j,i)
    132132                ENDIF
    133              ELSEIF ( icloud_scheme == 0  .AND.  .NOT. precipitation )  THEN
     133             ELSE
    134134                IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
    135135                   qc(k,j,i) = q(k,j,i) - q_s
    136136                   ql(k,j,i) = qc(k,j,i)
    137137                ELSE
    138                    qc(k,j,i) = 0.0_wp 
     138                   qc(k,j,i) = 0.0_wp
    139139                   ql(k,j,i) = 0.0_wp
    140                 ENDIF
    141              ELSE
    142                 IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
    143                    ql(k,j,i) = q(k,j,i) - q_s
    144                 ELSE
    145                    ql(k,j,i) = 0.0_wp
    146140                ENDIF
    147141             ENDIF
  • palm/trunk/SOURCE/check_parameters.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! PALM collision kernel deleted. Collision algorithms are checked for correct
     22! spelling.
     23!
     24! Tails removed.
     25!
     26! Checks for use_sgs_for_particles adopted for the use of droplets with
     27! use_sgs_for_particles adopted.
     28!
     29! Unused variables removed.
    2230!
    2331! Former revisions:
     
    358366    CHARACTER (LEN=1)   ::  sq                       !<
    359367    CHARACTER (LEN=15)  ::  var                      !<
    360     CHARACTER (LEN=7)   ::  unit                     !< 
     368    CHARACTER (LEN=7)   ::  unit                     !<
    361369    CHARACTER (LEN=8)   ::  date                     !<
    362370    CHARACTER (LEN=10)  ::  time                     !<
     
    367375
    368376    INTEGER(iwp) ::  i                               !<
    369     INTEGER(iwp) ::  ilen                            !<
    370     INTEGER(iwp) ::  iremote = 0                     !<
     377    INTEGER(iwp) ::  ilen                            !<
    371378    INTEGER(iwp) ::  j                               !<
    372379    INTEGER(iwp) ::  k                               !<
     
    374381    INTEGER(iwp) ::  netcdf_data_format_save         !<
    375382    INTEGER(iwp) ::  position                        !<
    376     INTEGER(iwp) ::  prec                            !<
    377383   
    378384    LOGICAL     ::  found                            !<
     
    661667                        'with particle advection.'
    662668       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
    663     ENDIF
    664 
    665 !
    666 !--
    667     IF ( use_particle_tails )  THEN
    668        message_string = 'Particle tails are currently not available due ' //   &
    669                         'to the new particle structure.'
    670        CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
    671669    ENDIF
    672670
     
    767765!
    768766!-- Check cloud scheme
    769     IF ( cloud_scheme == 'seifert_beheng' )  THEN
    770        icloud_scheme = 0
     767    IF ( cloud_scheme == 'saturation_adjust' )  THEN
     768       microphysics_sat_adjust = .TRUE.
     769       microphysics_seifert    = .FALSE.
     770       microphysics_kessler    = .FALSE.
     771       precipitation           = .FALSE.
     772    ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
     773       microphysics_sat_adjust = .FALSE.
     774       microphysics_seifert    = .TRUE.
     775       microphysics_kessler    = .FALSE.
     776       precipitation           = .TRUE.
    771777    ELSEIF ( cloud_scheme == 'kessler' )  THEN
    772        icloud_scheme = 1
     778       microphysics_sat_adjust = .FALSE.
     779       microphysics_seifert    = .FALSE.
     780       microphysics_kessler    = .TRUE.
     781       precipitation           = .TRUE.
    773782    ELSE
    774783       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
     
    846855    ENDIF
    847856
    848     IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke .AND.        &
    849          ( scalar_advec /= 'ws-scheme' .OR.                                    &
    850            scalar_advec /= 'ws-scheme-mono' )                                  &
     857    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.               &
     858         .NOT. use_upstream_for_tke  .AND.                                       &
     859         ( scalar_advec /= 'ws-scheme'  .OR.  scalar_advec /= 'ws-scheme-mono' ) &
    851860       )  THEN
    852861       use_upstream_for_tke = .TRUE.
     
    906915          hall_kernel = .TRUE.
    907916
    908        CASE ( 'palm' )
    909           palm_kernel = .TRUE.
    910 
    911917       CASE ( 'wang', 'wang_fast' )
    912918          wang_kernel = .TRUE.
     
    922928    END SELECT
    923929    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
     930
     931!
     932!-- Collision algorithms:
     933    SELECT CASE ( TRIM( collision_algorithm ) )
     934
     935       CASE ( 'all_or_nothing' )
     936          all_or_nothing = .TRUE.
     937
     938       CASE ( 'average_impact' )
     939          average_impact = .TRUE.
     940
     941       CASE DEFAULT
     942          message_string = 'unknown collision algorithm: collision_algorithm = "' // &
     943                           TRIM( collision_algorithm ) // '"'
     944          CALL message( 'check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
     945
     946       END SELECT
    924947
    925948    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     
    9791002              'not allowed with humidity = ', humidity
    9801003       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
    981     ENDIF
    982 
    983     IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
    984        WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ',   &
    985               'not allowed with cloud_physics = ', cloud_physics
    986        CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
    9871004    ENDIF
    9881005
     
    10321049       ENDIF
    10331050
    1034        IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
     1051       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    10351052          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
    10361053                           ' seifert_beheng'
     
    10951112    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
    10961113                 loop_optimization == 'vector' )                               &
    1097          .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
     1114         .AND.  cloud_physics  .AND.  microphysics_seifert )  THEN
    10981115       message_string = 'cloud_scheme = seifert_beheng requires ' //           &
    10991116                        'loop_optimization = "cache" or "vector"'
     
    21522169!
    21532170!-- Set the default value for the integration interval of precipitation amount
    2154     IF ( precipitation )  THEN
     2171    IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
    21552172       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
    21562173          precipitation_amount_interval = dt_do2d_xy
     
    27822799                                 'lemented for cloud_physics = .FALSE.'
    27832800                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
    2784              ELSEIF ( icloud_scheme /= 0 )  THEN
     2801             ELSEIF ( .NOT.  microphysics_seifert )  THEN
    27852802                message_string = 'data_output_pr = ' //                        &
    27862803                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    27872804                                 'lemented for cloud_scheme /= seifert_beheng'
    27882805                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    2789              ELSEIF (  .NOT.  precipitation )  THEN
    2790                 message_string = 'data_output_pr = ' //                        &
    2791                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
    2792                                  'lemented for precipitation = .FALSE.'
    2793                 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
    27942806             ELSE
    27952807                dopr_index(i) = 73
     
    28042816                                 'lemented for cloud_physics = .FALSE.'
    28052817                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
    2806              ELSEIF ( icloud_scheme /= 0 )  THEN
     2818             ELSEIF ( .NOT.  microphysics_seifert )  THEN
    28072819                message_string = 'data_output_pr = ' //                        &
    28082820                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    28092821                                 'lemented for cloud_scheme /= seifert_beheng'
    28102822                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    2811              ELSEIF (  .NOT.  precipitation )  THEN
    2812                 message_string = 'data_output_pr = ' //                        &
    2813                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
    2814                                  'lemented for precipitation = .FALSE.'
    2815                 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
    28162823             ELSE
    28172824                dopr_index(i) = 74
     
    28262833                                 'lemented for cloud_physics = .FALSE.'
    28272834                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
    2828              ELSEIF ( icloud_scheme /= 0 )  THEN
    2829                 message_string = 'data_output_pr = ' //                        &
    2830                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
    2831                                  'lemented for cloud_scheme /= seifert_beheng'
    2832                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    28332835             ELSE
    28342836                dopr_index(i) = 75
     
    28432845                                 'lemented for cloud_physics = .FALSE.'
    28442846                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
    2845              ELSEIF ( icloud_scheme /= 0 )  THEN
    2846                 message_string = 'data_output_pr = ' //                        &
    2847                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
    2848                                  'lemented for cloud_scheme /= seifert_beheng'
    2849                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    2850              ELSEIF (  .NOT.  precipitation )  THEN
    2851                 message_string = 'data_output_pr = ' //                        &
    2852                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
    2853                                  'lemented for precipitation = .FALSE.'
    2854                 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
    2855 
     2847             ELSEIF ( microphysics_sat_adjust )  THEN
     2848                message_string = 'data_output_pr = ' //                        &
     2849                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     2850                                 'ilable for cloud_scheme = saturation_adjust'
     2851                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
    28562852             ELSE
    28572853                dopr_index(i) = 76
     
    32033199                         'res cloud_physics = .TRUE.'
    32043200                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    3205              ELSEIF ( icloud_scheme /= 0 )  THEN
     3201             ELSEIF ( .NOT.  microphysics_seifert )  THEN
    32063202                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    32073203                         'res cloud_scheme = seifert_beheng'
     
    32243220                         'res cloud_physics = .TRUE.'
    32253221                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    3226              ELSEIF ( icloud_scheme /= 0 )  THEN
    3227                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3228                          'res cloud_scheme = seifert_beheng'
    3229                 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    3230              ELSEIF (  .NOT.  precipitation )  THEN
    3231                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3232                                  'res precipitation = .TRUE.'
    3233                 CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
     3222             ELSEIF ( microphysics_sat_adjust )  THEN
     3223                message_string = 'output of "' // TRIM( var ) // '" is ' //  &
     3224                         'not available for cloud_scheme = saturation_adjust'
     3225                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
    32343226             ENDIF
    32353227             unit = 'kg/kg m/s'
     
    32493241                         'res cloud_physics = .TRUE.'
    32503242                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    3251              ELSEIF ( icloud_scheme /= 0 ) THEN
    3252                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3253                          'res cloud_scheme = seifert_beheng'
    3254                 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    32553243             ENDIF
    32563244             unit = 'kg/kg'
     
    32793267                         'res cloud_physics = .TRUE.'
    32803268                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    3281              ELSEIF ( icloud_scheme /= 0 ) THEN
     3269             ELSEIF ( .NOT.  microphysics_seifert ) THEN
    32823270                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    32833271                         'res cloud_scheme = seifert_beheng'
    32843272                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    3285              ELSEIF (  .NOT.  precipitation )  THEN
    3286                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3287                                  'res precipitation = .TRUE.'
    3288                 CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    32893273             ENDIF
    32903274             unit = 'kg/kg'
     
    33713355                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    33723356             ENDIF
    3373              IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
     3357             IF ( TRIM( var ) == 'pra*'  .AND.                                 &
     3358                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
    33743359                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3375                                  'res precipitation = .TRUE.'
     3360                                 'res cloud_scheme = kessler or seifert_beheng'
    33763361                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    33773362             ENDIF
     
    33813366                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
    33823367             ENDIF
    3383              IF ( TRIM( var ) == 'prr*'  .AND.  .NOT.  precipitation )  THEN
     3368             IF ( TRIM( var ) == 'prr*'  .AND.                                 &
     3369                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
    33843370                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3385                                  'res precipitation = .TRUE.'
     3371                                 'res cloud_scheme = kessler or seifert_beheng'
    33863372                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    33873373             ENDIF
  • palm/trunk/SOURCE/data_output_2d.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Output of bulk cloud physics simplified.
    2222!
    2323! Former revisions:
     
    146146               do2d_xz_last_time, do2d_xz_n, do2d_xz_time_count,               &
    147147               do2d_yz_last_time, do2d_yz_n, do2d_yz_time_count,               &
    148                ibc_uv_b, icloud_scheme, io_blocks, io_group,                   &
    149                message_string,                             &
    150                ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz, psolver, section,        &
    151                simulated_time,  simulated_time_chr, time_since_reference_point
     148               ibc_uv_b, io_blocks, io_group, message_string,                  &
     149               ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
     150               psolver, section, simulated_time, simulated_time_chr,           &
     151               time_since_reference_point
    152152       
    153153    USE cpulog,                                                                &
     
    683683
    684684             CASE ( 'prr*_xy' )        ! 2d-array
    685                 IF ( icloud_scheme == 1 )  THEN
    686                    IF ( av == 0 )  THEN
    687                       CALL exchange_horiz_2d( precipitation_rate )
    688                       DO  i = nxlg, nxrg
    689                          DO  j = nysg, nyng
    690                             local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
     685                IF ( av == 0 )  THEN
     686                   CALL exchange_horiz_2d( prr(nzb+1,:,:) )
     687                   DO  i = nxlg, nxrg
     688                      DO  j = nysg, nyng
     689                         local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
     690                      ENDDO
     691                   ENDDO
     692                ELSE
     693                   CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
     694                   DO  i = nxlg, nxrg
     695                      DO  j = nysg, nyng
     696                         local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) * hyrho(nzb+1)
     697                      ENDDO
     698                   ENDDO
     699                ENDIF
     700                resorted = .TRUE.
     701                two_d = .TRUE.
     702                level_z(nzb+1) = zu(nzb+1)
     703
     704             CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
     705                IF ( av == 0 )  THEN
     706                   CALL exchange_horiz( prr, nbgp )
     707                   DO  i = nxlg, nxrg
     708                      DO  j = nysg, nyng
     709                         DO  k = nzb, nzt+1
     710                            local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
    691711                         ENDDO
    692712                      ENDDO
    693                    ELSE
    694                       CALL exchange_horiz_2d( precipitation_rate_av )
    695                       DO  i = nxlg, nxrg
    696                          DO  j = nysg, nyng
    697                             local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
    698                          ENDDO
    699                       ENDDO
    700                    ENDIF
    701                 ELSE
    702                    IF ( av == 0 )  THEN
    703                       CALL exchange_horiz_2d( prr(nzb+1,:,:) )
    704                       DO  i = nxlg, nxrg
    705                          DO  j = nysg, nyng
    706                             local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
    707                          ENDDO
    708                       ENDDO
    709                    ELSE
    710                       CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
    711                       DO  i = nxlg, nxrg
    712                          DO  j = nysg, nyng
    713                             local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) *          &
    714                                                   hyrho(nzb+1)
    715                          ENDDO
    716                       ENDDO
    717                    ENDIF
    718                 ENDIF
    719                 resorted = .TRUE.
    720                 two_d = .TRUE.
    721                 level_z(nzb+1) = zu(nzb+1)
    722 
    723              CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
    724                 IF ( av == 0 )  THEN
    725                    CALL exchange_horiz( prr, nbgp )
     713                   ENDDO
     714                ELSE
     715                   CALL exchange_horiz( prr_av, nbgp )
    726716                   DO  i = nxlg, nxrg
    727717                      DO  j = nysg, nyng
    728718                         DO  k = nzb, nzt+1
    729                             local_pf(i,j,k) = prr(k,j,i)
    730                          ENDDO
    731                       ENDDO
    732                    ENDDO
    733                 ELSE
    734                    CALL exchange_horiz( prr_av, nbgp )
    735                    DO  i = nxlg, nxrg
    736                       DO  j = nysg, nyng
    737                          DO  k = nzb, nzt+1
    738                             local_pf(i,j,k) = prr_av(k,j,i)
     719                            local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
    739720                         ENDDO
    740721                      ENDDO
  • palm/trunk/SOURCE/data_output_3d.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! prr vertical dimensions set to nzb_do to nzt_do. Unused variables deleted.
    2222!
    2323! Former revisions:
     
    122122       
    123123    USE control_parameters,                                                    &
    124         ONLY:  avs_data_file, cloud_physics, do3d, do3d_avs_n,                 &
    125                do3d_no, do3d_time_count, io_blocks, io_group,                  &
    126                message_string, ntdim_3d,                                       &
    127                nz_do3d, plot_3d_precision, psolver, simulated_time,            &
    128                simulated_time_chr, skip_do_avs, time_since_reference_point
     124        ONLY:  cloud_physics, do3d, do3d_no, do3d_time_count, io_blocks,       &
     125               io_group, message_string, ntdim_3d, nz_do3d, psolver,           &
     126               simulated_time, time_since_reference_point
    129127       
    130128    USE cpulog,                                                                &
     
    132130       
    133131    USE indices,                                                               &
    134         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzt,  &
    135                nzb
     132        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb
    136133       
    137134    USE kinds
     
    163160    IMPLICIT NONE
    164161
    165     CHARACTER (LEN=9) ::  simulated_time_mod  !<
    166 
    167162    INTEGER(iwp) ::  av        !<
    168163    INTEGER(iwp) ::  i         !<
     
    173168    INTEGER(iwp) ::  nzb_do    !< vertical lower limit for data output
    174169    INTEGER(iwp) ::  nzt_do    !< vertical upper limit for data output
    175     INTEGER(iwp) ::  pos       !<
    176     INTEGER(iwp) ::  prec      !<
    177     INTEGER(iwp) ::  psi       !<
    178170
    179171    LOGICAL      ::  found     !<
     
    374366                DO  i = nxlg, nxrg
    375367                   DO  j = nysg, nyng
    376                       DO  k = nzb, nzt+1
     368                      DO  k = nzb_do, nzt_do
    377369                         local_pf(i,j,k) = prr(k,j,i)
    378370                      ENDDO
     
    383375                DO  i = nxlg, nxrg
    384376                   DO  j = nysg, nyng
    385                       DO  k = nzb, nzt+1
     377                      DO  k = nzb_do, nzt_do
    386378                         local_pf(i,j,k) = prr_av(k,j,i)
    387379                      ENDDO
  • palm/trunk/SOURCE/data_output_dvrp.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Particles and tails removed.
    2222!
    2323! Former revisions:
     
    5858! Description:
    5959! ------------
    60 !> Plot of isosurface, particles and slicers with dvrp-software
     60!> Plot of isosurface and slicers with dvrp-software
    6161!------------------------------------------------------------------------------!
    6262 MODULE dvrp_color
     
    130130       
    131131    USE kinds
    132    
    133     USE particle_attributes,                                                   &
    134         ONLY:  maximum_number_of_tailpoints, number_of_particles,              &
    135                number_of_tails, particle_advection, particle_advection_start,  &
    136                particle_tail_coordinates, particles, uniform_particles,        &
    137                use_particle_tails
    138        
     132
    139133    USE pegrid
    140134
     
    234228!--    Select the plot mode (in case of isosurface or slicer only if user has
    235229!--    defined a variable which shall be plotted; otherwise do nothing)
    236        IF ( mode_dvrp(m)(1:9) == 'particles'  .AND.  particle_advection  .AND. &
    237             simulated_time >= particle_advection_start )  THEN
    238 !
    239 !--       DVRP-Calls for plotting particles:
    240           CALL cpu_log( log_point_s(28), 'dvrp_particles', 'start' )
    241 
    242 !
    243 !--       Definition of characteristics of particle material
    244 !          tmp_r = 0.1;  tmp_g = 0.7;  tmp_b = 0.1;  tmp_t = 0.0
    245           tmp_r = 0.0_wp;  tmp_g = 0.0_wp;  tmp_b = 0.0_wp;  tmp_t = 0.0_wp
    246           CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
    247 
    248 !
    249 !--       If clipping is active and if this subdomain is clipped, find out the
    250 !--       number of particles and tails to be plotted; otherwise, all
    251 !--       particles/tails are plotted. dvrp_mask is used to mark the partikles.
    252           ALLOCATE( dvrp_mask(number_of_particles) )
    253 
    254           IF ( dvrp_total_overlap )  THEN
    255              dvrp_mask = .TRUE.
    256              dvrp_nop  = number_of_particles
    257              dvrp_not  = number_of_tails
    258           ELSE
    259              dvrp_mask = .FALSE.
    260              dvrp_nop  = 0
    261              dvrp_not  = 0
    262              IF ( dvrp_overlap )  THEN
    263                 IF ( .NOT. use_particle_tails )  THEN
    264                    DO  n = 1, number_of_particles
    265                       ip = particles(n)%x / dx
    266                       jp = particles(n)%y / dy
    267                       IF ( ip >= nxl_dvrp  .AND.  ip <= nxr_dvrp  .AND.  &
    268                            jp >= nys_dvrp  .AND.  jp <= nyn_dvrp )  THEN
    269                          dvrp_nop = dvrp_nop + 1
    270                          dvrp_mask(n) = .TRUE.
    271                       ENDIF
    272                    ENDDO
    273                 ELSE
    274                    k = 0
    275                    DO  n = 1, number_of_particles
    276                       IF ( particles(n)%tail_id /= 0 )  THEN
    277                          k = k + 1
    278                          ip = particles(n)%x / dx
    279                          jp = particles(n)%y / dy
    280                          IF ( ip >= nxl_dvrp  .AND.  ip <= nxr_dvrp  .AND.  &
    281                               jp >= nys_dvrp  .AND.  jp <= nyn_dvrp )  THEN
    282                             dvrp_not = dvrp_not + 1
    283                             dvrp_mask(n) = .TRUE.
    284                          ENDIF
    285                       ENDIF
    286                    ENDDO
    287                 ENDIF
    288              ENDIF
    289           ENDIF
    290 
    291 !
    292 !--       Move particle coordinates to one-dimensional arrays
    293           IF ( .NOT. use_particle_tails )  THEN
    294 !
    295 !--          All particles are output
    296              ALLOCATE( psize(dvrp_nop), p_t(dvrp_nop), p_c(dvrp_nop), &
    297                        p_x(dvrp_nop), p_y(dvrp_nop), p_z(dvrp_nop) )
    298              psize = 0.0_wp;  p_t = 0_wp;  p_c = 0.0_wp
    299              p_x = 0.0_wp; p_y = 0.0_wp
    300              p_z   = 0.0_wp
    301              k = 0
    302              DO  n = 1, number_of_particles
    303                 IF ( dvrp_mask(n) )  THEN
    304                    k = k + 1
    305                    psize(k) = particles(n)%dvrp_psize
    306                    p_x(k)   = particles(n)%x * superelevation_x
    307                    p_y(k)   = particles(n)%y * superelevation_y
    308                    p_z(k)   = particles(n)%z * superelevation
    309                    p_c(k)   = particles(n)%class
    310                 ENDIF
    311              ENDDO
    312           ELSE
    313 !
    314 !--          Particles have a tail
    315              ALLOCATE( psize(dvrp_not), p_t(dvrp_not),             &
    316                        p_c(dvrp_not*maximum_number_of_tailpoints), &
    317                        p_x(dvrp_not*maximum_number_of_tailpoints), &
    318                        p_y(dvrp_not*maximum_number_of_tailpoints), &
    319                        p_z(dvrp_not*maximum_number_of_tailpoints) )
    320              psize = 0.0_wp;  p_t = 0_wp;  p_c = 0.0_wp
    321              p_x = 0.0_wp;  p_y = 0.0_wp
    322              p_z   = 0.0_wp
    323              i = 0
    324              k = 0
    325 
    326              DO  n = 1, number_of_particles
    327                 nn = particles(n)%tail_id
    328                 IF ( nn /= 0  .AND.  dvrp_mask(n) )  THEN
    329                    k = k + 1
    330                    DO  j = 1, particles(n)%tailpoints
    331                       i = i + 1
    332                       p_x(i) = particle_tail_coordinates(j,1,nn) * &
    333                                                                 superelevation_x
    334                       p_y(i) = particle_tail_coordinates(j,2,nn) * &
    335                                                                 superelevation_y
    336                       p_z(i) = particle_tail_coordinates(j,3,nn) * &
    337                                                                 superelevation
    338                       p_c(i) = particle_tail_coordinates(j,4,nn)
    339                    ENDDO
    340                    psize(k) = particles(n)%dvrp_psize
    341                    p_t(k)   = particles(n)%tailpoints - 1
    342                 ENDIF               
    343              ENDDO
    344           ENDIF
    345 
    346 !
    347 !--       Compute and plot particles in dvr-format
    348           IF ( uniform_particles  .AND.  .NOT. use_particle_tails )  THEN
    349 !
    350 !--          All particles have the same color. Use simple routine to set
    351 !--          the particle attributes (produces less output data)
    352              CALL DVRP_PARTICLES( m-1, p_x, p_y, p_z, psize )
    353           ELSE
    354 !
    355 !--          Set color definitions
    356              CALL user_dvrp_coltab( 'particles', 'none' )
    357 
    358              CALL DVRP_COLORTABLE_HLS( m-1, 0, interval_values_dvrp_prt, &
    359                                        interval_h_dvrp_prt,              &
    360                                        interval_l_dvrp_prt,              &
    361                                        interval_s_dvrp_prt,              &
    362                                        interval_a_dvrp_prt )
    363 
    364              IF ( .NOT. use_particle_tails )  THEN
    365                 CALL DVRP_PARTICLES( m-1, dvrp_nop, p_x, p_y, p_z, 3, psize, &
    366                                      p_c, p_t )
    367              ELSE
    368                 CALL DVRP_PARTICLES( m-1, dvrp_not, p_x, p_y, p_z, 15, psize, &
    369                                      p_c, p_t )
    370              ENDIF
    371           ENDIF
    372 
    373           CALL DVRP_VISUALIZE( m-1, 3, dvrp_filecount )
    374 
    375           DEALLOCATE( dvrp_mask, psize, p_c, p_t, p_x, p_y, p_z )
    376 
    377           CALL cpu_log( log_point_s(28), 'dvrp_particles', 'stop' )
    378 
    379 
    380        ELSEIF ( ( mode_dvrp(m)(1:10) == 'isosurface'  .OR.                     &
     230       IF ( ( mode_dvrp(m)(1:10) == 'isosurface'  .OR.                     &
    381231                  mode_dvrp(m)(1:6)  == 'slicer'           )                   &
    382232                  .AND.  output_variable /= ' ' )  THEN
  • palm/trunk/SOURCE/flow_statistics.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Output of bulk microphysics simplified.
    2222!
    2323! Former revisions:
     
    193193    USE control_parameters,                                                    &
    194194        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    195                 dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
     195                dt_3d, g, humidity, kappa, large_scale_forcing,                &
    196196                large_scale_subsidence, max_pr_user, message_string, neutral,  &
    197                 ocean, passive_scalar, precipitation, simulated_time,          &
     197                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
    198198                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    199199                ws_scheme_mom, ws_scheme_sca
     
    898898                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
    899899                                                          rmask(j,i,sr)
     900                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
     901
    900902                      IF ( .NOT. cloud_droplets )  THEN
    901903                         pts = 0.5_wp *                                        &
     
    906908                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
    907909                                                             rmask(j,i,sr)
    908                          IF ( icloud_scheme == 0  )  THEN
    909                             sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
     910                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
     911                                                             rmask(j,i,sr)
     912                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
     913                                                             rmask(j,i,sr)
     914                         IF ( microphysics_seifert )  THEN
     915                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
    910916                                                                rmask(j,i,sr)
    911                             sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *    &
    912                                                                 rmask(j,i,sr)
    913                             IF ( precipitation )  THEN
    914                                sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
    915                                                                    rmask(j,i,sr)
    916                                sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
    917                                                                    rmask(j,i,sr)
    918                                sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
    919                                                                    rmask(j,i,sr)
    920                             ENDIF
    921                          ELSE
    922                             sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
     917                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
    923918                                                                rmask(j,i,sr)
    924919                         ENDIF
    925                       ELSE
    926                          sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *       &
    927                                                              rmask(j,i,sr)
    928920                      ENDIF
     921
    929922                   ELSE
    930923                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
     
    16621655    USE control_parameters,                                                    &
    16631656        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    1664                 dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
    1665                 large_scale_subsidence, max_pr_user, message_string, neutral,  &
    1666                 ocean, passive_scalar, precipitation, simulated_time,          &
    1667                 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    1668                 ws_scheme_mom, ws_scheme_sca
     1657                dt_3d, g, humidity, kappa, large_scale_forcing,                &
     1658                large_scale_subsidence, max_pr_user, message_string,           &
     1659                microphysics_seifert, neutral, ocean, passive_scalar,          &
     1660                simulated_time, use_subsidence_tendencies, use_surface_fluxes, &
     1661                use_top_fluxes, ws_scheme_mom, ws_scheme_sca
    16691662       
    16701663    USE cpulog,                                                                &
     
    27742767                !$acc end parallel loop
    27752768
    2776                 IF ( icloud_scheme == 0 )  THEN
     2769                IF ( microphysics_seifert )  THEN
    27772770
    27782771                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
     
    27912784                   !$acc end parallel loop
    27922785
    2793                    IF ( precipitation )  THEN
    2794 
    2795                       !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
    2796                       DO  k = nzb, nzt_diff
    2797                          s1 = 0
    2798                          s2 = 0
    2799                          s3 = 0
    2800                          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    2801                          DO  i = nxl, nxr
    2802                             DO  j = nys, nyn
    2803                                s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2804                                s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2805                                s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2806                             ENDDO
     2786                   !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
     2787                   DO  k = nzb, nzt_diff
     2788                      s1 = 0
     2789                      s2 = 0
     2790                      s3 = 0
     2791                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
     2792                      DO  i = nxl, nxr
     2793                         DO  j = nys, nyn
     2794                            s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2795                            s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2796                            s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    28072797                         ENDDO
    2808                          sums_l(k,73,tn) = s1
    2809                          sums_l(k,74,tn) = s2
    2810                          sums_l(k,76,tn) = s3
    28112798                      ENDDO
    2812                       !$acc end parallel loop
    2813 
    2814                    ENDIF
     2799                      sums_l(k,73,tn) = s1
     2800                      sums_l(k,74,tn) = s2
     2801                      sums_l(k,76,tn) = s3
     2802                   ENDDO
     2803                   !$acc end parallel loop
    28152804
    28162805                ELSE
  • palm/trunk/SOURCE/header.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
    22 ! 
     21! Tails removed. icloud_scheme replaced by microphysics_*
     22!
    2323! Former revisions:
    2424! -----------------
     
    275275               density_ratio, dissipation_classes, dt_min_part, dt_prel,       &
    276276               dt_write_particle_data, end_time_prel,                          &
    277                maximum_number_of_tailpoints, maximum_tailpoint_age,            &
    278                minimum_tailpoint_distance, number_of_particle_groups,          &
    279                particle_advection, particle_advection_start,                   &
     277               number_of_particle_groups, particle_advection,                  &
     278               particle_advection_start,                                       &
    280279               particles_per_point, pdx, pdy, pdz,  psb, psl, psn, psr, pss,   &
    281280               pst, radius, radius_classes, random_start_position,             &
    282281               seed_follows_topography,                                        &
    283                total_number_of_particles, use_particle_tails,                  &
    284                use_sgs_for_particles, total_number_of_tails,                   &
     282               total_number_of_particles, use_sgs_for_particles,               &
    285283               vertical_particle_advection, write_particle_statistics
    286284       
     
    17001698                                   slicer_range_limits_dvrp(:,m)
    17011699             ENDIF
    1702           ELSEIF ( mode_dvrp(i)(1:9) == 'particles' )  THEN
    1703              WRITE ( io, 363 )  dvrp_psize
    1704              IF ( particle_dvrpsize /= 'none' )  THEN
    1705                 WRITE ( io, 364 )  'size', TRIM( particle_dvrpsize ), &
    1706                                    dvrpsize_interval
    1707              ENDIF
    1708              IF ( particle_color /= 'none' )  THEN
    1709                 WRITE ( io, 364 )  'color', TRIM( particle_color ), &
    1710                                    color_interval
    1711              ENDIF
    17121700          ENDIF
    17131701          i = i + 1
     
    17821770       WRITE ( io, 415 )
    17831771       WRITE ( io, 416 ) surface_pressure, r_d, rho_surface, cp, l_v
    1784        IF ( icloud_scheme == 0 )  THEN
     1772       IF ( microphysics_seifert )  THEN
    17851773          WRITE ( io, 510 ) 1.0E-6_wp * nc_const
    1786           IF ( precipitation )  WRITE ( io, 511 ) c_sedimentation
     1774          WRITE ( io, 511 ) c_sedimentation
    17871775       ENDIF
    17881776    ENDIF
     
    17961784       WRITE ( io, 432 )
    17971785       IF ( cloud_top_radiation )  WRITE ( io, 132 )
    1798        IF ( icloud_scheme == 1 )  THEN
    1799           IF ( precipitation )  WRITE ( io, 133 )
    1800        ELSEIF ( icloud_scheme == 0 )  THEN
     1786       IF ( microphysics_kessler )  THEN
     1787          WRITE ( io, 133 )
     1788       ELSEIF ( microphysics_seifert )  THEN
    18011789          IF ( drizzle )  WRITE ( io, 506 )
    1802           IF ( precipitation )  THEN
    1803              WRITE ( io, 505 )
    1804              IF ( turbulence )  WRITE ( io, 507 )
    1805              IF ( ventilation_effect )  WRITE ( io, 508 )
    1806              IF ( limiter_sedimentation )  WRITE ( io, 509 )
    1807           ENDIF
     1790          WRITE ( io, 505 )
     1791          IF ( turbulence )  WRITE ( io, 507 )
     1792          IF ( ventilation_effect )  WRITE ( io, 508 )
     1793          IF ( limiter_sedimentation )  WRITE ( io, 509 )
    18081794       ENDIF
    18091795    ELSEIF ( humidity  .AND.  cloud_droplets )  THEN
     
    18721858       IF ( particles_per_point > 1 )  WRITE ( io, 489 )  particles_per_point
    18731859       WRITE ( io, 495 )  total_number_of_particles
    1874        IF ( use_particle_tails  .AND.  maximum_number_of_tailpoints /= 0 )  THEN
    1875           WRITE ( io, 483 )  maximum_number_of_tailpoints
    1876           IF ( minimum_tailpoint_distance /= 0 )  THEN
    1877              WRITE ( io, 484 )  total_number_of_tails,      &
    1878                                 minimum_tailpoint_distance, &
    1879                                 maximum_tailpoint_age
    1880           ENDIF
    1881        ENDIF
    18821860       IF ( dt_write_particle_data /= 9999999.9_wp )  THEN
    18831861          WRITE ( io, 485 )  dt_write_particle_data
     
    22152193362 FORMAT (/'       Slicer plane ',A/ &
    22162194            '       Slicer limits: [',F6.2,',',F6.2,']')
    2217 363 FORMAT (/'       Particles'/ &
    2218             '          particle size:  ',F7.2,' m')
    2219 364 FORMAT ('          particle ',A,' controlled by "',A,'" with interval [', &
    2220                        F6.2,',',F6.2,']')
    22212195365 FORMAT (/'       Groundplate color: (',F4.2,',',F4.2,',',F4.2,') (R,G,B)'/ &
    22222196            '       Superelevation along (x,y,z): (',F4.1,',',F4.1,',',F4.1, &
     
    23772351481 FORMAT ('       Particles have random start positions'/)
    23782352482 FORMAT ('          Particles are advected only horizontally'/)
    2379 483 FORMAT ('       Particles have tails with a maximum of ',I3,' points')
    2380 484 FORMAT ('            Number of tails of the total domain: ',I10/ &
    2381             '            Minimum distance between tailpoints: ',F8.2,' m'/ &
    2382             '            Maximum age of the end of the tail:  ',F8.2,' s')
    23832353485 FORMAT ('       Particle data are written on file every ', F9.1, ' s')
    23842354486 FORMAT ('       Particle statistics are written on file'/)
  • palm/trunk/SOURCE/init_3d_model.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! icloud_scheme replaced by microphysics_*
    2222!
    2323! Former revisions:
     
    465465                       precipitation_rate(nysg:nyng,nxlg:nxrg) )
    466466
    467              IF ( icloud_scheme == 0 )  THEN
    468 !
    469 !--             1D-arrays
    470                 ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1),                 &
    471                            q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) )
    472 !
    473 !--             3D-cloud water content
     467!
     468!--          1D-arrays
     469             ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1),                 &
     470                        q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) )
     471!
     472!--          3D-cloud water content
    474473#if defined( __nopointer )
    475                 ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     474             ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    476475#else
    477                 ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     476             ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    478477#endif
    479 
    480                 IF ( precipitation )  THEN
    481 !
    482 !--                1D-arrays
    483                    ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
    484 
    485 !
    486 !--                2D-rain water content and rain drop concentration arrays
    487                    ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
    488                               qrsws(nysg:nyng,nxlg:nxrg),                      &
    489                               qrswst(nysg:nyng,nxlg:nxrg),                     &
    490                               nrs(nysg:nyng,nxlg:nxrg),                        &
    491                               nrsws(nysg:nyng,nxlg:nxrg),                      &
    492                               nrswst(nysg:nyng,nxlg:nxrg) )
    493 !
    494 !--                3D-rain water content, rain drop concentration arrays
     478!
     479!--          3d-precipitation rate
     480             ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     481
     482             IF ( microphysics_seifert )  THEN
     483!
     484!--             1D-arrays
     485                ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
     486
     487!
     488!--             2D-rain water content and rain drop concentration arrays
     489                ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
     490                           qrsws(nysg:nyng,nxlg:nxrg),                      &
     491                           qrswst(nysg:nyng,nxlg:nxrg),                     &
     492                           nrs(nysg:nyng,nxlg:nxrg),                        &
     493                           nrsws(nysg:nyng,nxlg:nxrg),                      &
     494                           nrswst(nysg:nyng,nxlg:nxrg) )
     495!
     496!--             3D-rain water content, rain drop concentration arrays
    495497#if defined( __nopointer )
    496                    ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
    497                              nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    498                              qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
    499                              qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    500                              tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
    501                              tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     498                ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
     499                          nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     500                          qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
     501                          qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     502                          tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
     503                          tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    502504#else
    503                    ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    504                              nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    505                              nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    506                              qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    507                              qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    508                              qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     505                ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     506                          nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     507                          nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     508                          qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     509                          qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     510                          qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    509511#endif
    510 !
    511 !--                3d-precipitation rate
    512                    ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    513                 ENDIF
    514 
    515              ENDIF
     512             ENDIF
     513
    516514          ENDIF
    517515
     
    658656          IF ( cloud_physics )  THEN
    659657             ql => ql_1
    660              IF ( icloud_scheme == 0 )  THEN
    661                 qc => qc_1
    662                 IF ( precipitation )  THEN
    663                    qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
    664                    nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
    665                 ENDIF
     658             qc => qc_1
     659             IF ( microphysics_seifert )  THEN
     660                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
     661                nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
    666662             ENDIF
    667663          ENDIF
     
    733729                ENDDO
    734730             ENDDO
    735              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    736                   precipitation )  THEN
     731             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    737732                DO  i = nxlg, nxrg
    738733                   DO  j = nysg, nyng
     
    790785          IF ( humidity  .OR.  passive_scalar )  THEN
    791786             qs = 0.0_wp
    792              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    793                   precipitation )  THEN
     787             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    794788                qrs = 0.0_wp
    795789                nrs = 0.0_wp
     
    880874                ENDDO
    881875             ENDDO
    882              IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
     876             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    883877!
    884878!--             Initialze nc_1d with default value
    885879                nc_1d(:) = nc_const
    886880
    887                 IF ( precipitation )  THEN
    888                    DO  i = nxlg, nxrg
    889                       DO  j = nysg, nyng
    890                          qr(:,j,i) = 0.0_wp
    891                          nr(:,j,i) = 0.0_wp
    892                       ENDDO
     881                DO  i = nxlg, nxrg
     882                   DO  j = nysg, nyng
     883                      qr(:,j,i) = 0.0_wp
     884                      nr(:,j,i) = 0.0_wp
    893885                   ENDDO
    894                 ENDIF
     886                ENDDO
    895887
    896888             ENDIF
     
    10781070!--       Determine the near-surface water flux
    10791071          IF ( humidity  .OR.  passive_scalar )  THEN
    1080              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    1081                   precipitation )  THEN
     1072             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    10821073                qrsws = 0.0_wp
    10831074                nrsws = 0.0_wp
     
    11161107             IF ( humidity  .OR.  passive_scalar )  THEN
    11171108                qswst = 0.0_wp
    1118                 IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.           &
    1119                      precipitation ) THEN
     1109                IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
    11201110                   nrswst = 0.0_wp
    11211111                   qrswst = 0.0_wp
     
    11571147          IF ( humidity  .OR.  passive_scalar )  THEN
    11581148             IF (  .NOT.  constant_waterflux )  qsws   = 0.0_wp
    1159              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    1160                   precipitation )  THEN
     1149             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    11611150                qrsws = 0.0_wp
    11621151                nrsws = 0.0_wp
     
    12001189       IF ( cloud_physics )  THEN
    12011190          ql = 0.0_wp
    1202           IF ( precipitation )  precipitation_amount = 0.0_wp
    1203           IF ( icloud_scheme == 0 )  THEN
    1204              qc = 0.0_wp
     1191          qc = 0.0_wp
     1192
     1193          precipitation_amount = 0.0_wp
     1194
     1195          IF ( microphysics_seifert )  THEN
    12051196             nc_1d = nc_const
    12061197          ENDIF
     
    12401231          tq_m = 0.0_wp
    12411232          q_p = q
    1242           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    1243                precipitation )  THEN
     1233          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    12441234             tqr_m = 0.0_wp
    1245              qr_p = qr
     1235             qr_p  = qr
    12461236             tnr_m = 0.0_wp
    1247              nr_p = nr
     1237             nr_p  = nr
    12481238          ENDIF
    12491239       ENDIF
     
    14261416       IF ( humidity  .OR.  passive_scalar )  THEN
    14271417          q_p = q
    1428           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    1429                precipitation )  THEN
     1418          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    14301419             qr_p = qr
    14311420             nr_p = nr
     
    14411430       IF ( humidity  .OR.  passive_scalar )  THEN
    14421431          tq_m = 0.0_wp
    1443           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    1444                precipitation )  THEN
     1432          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    14451433             tqr_m = 0.0_wp
    14461434             tnr_m = 0.0_wp
  • palm/trunk/SOURCE/init_cloud_physics.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! icloud_scheme replaced by microphysics_*
    2222!
    2323! Former revisions:
     
    9292       
    9393    USE control_parameters,                                                    &
    94         ONLY:  g, icloud_scheme, message_string, precipitation, pt_surface,    &
     94        ONLY:  g, message_string, microphysics_seifert, pt_surface,            &
    9595               rho_surface, surface_pressure
    9696   
     
    126126!
    127127!-- Calculate timestep according to precipitation
    128     IF ( icloud_scheme == 0  .AND.  precipitation )  THEN
     128    IF ( microphysics_seifert )  THEN
    129129       dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /         &
    130130                          w_precipitation
  • palm/trunk/SOURCE/init_masks.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! icloud_scheme replaced by microphysics_*
    2222!
    2323! Former revisions:
     
    9797               data_output_masks, data_output_masks_user,                      &
    9898               doav, doav_n, domask, domask_no, dz, dz_stretch_level, humidity,&
    99                icloud_scheme, mask, masks, mask_scale, mask_i,                 &
     99               mask, masks, mask_scale, mask_i,                                &
    100100               mask_i_global, mask_j, mask_j_global, mask_k, mask_k_global,    &
    101101               mask_loop, mask_size, mask_size_l, mask_start_l, mask_x,        &
    102102               mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
    103103               mask_z_loop, max_masks,  message_string, mid,                   &
    104                passive_scalar, precipitation, ocean
     104               microphysics_seifert, passive_scalar, ocean
    105105
    106106    USE grid_variables,                                                        &
     
    270270                        '" requires cloud_physics = .TRUE.'
    271271                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
    272                  ELSEIF ( icloud_scheme /= 0 ) THEN
     272                 ELSEIF ( .NOT. microphysics_seifert ) THEN
    273273                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    274274                         'res cloud_scheme = seifert_beheng'
    275275                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    276                  ELSEIF ( .NOT. precipitation )  THEN
    277                    message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    278                                  'res precipitation = .TRUE.'
    279                    CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    280276                ENDIF
    281277                unit = '1/m3'
     
    305301                            'res cloud_physics = .TRUE.'
    306302                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    307                 ELSEIF ( icloud_scheme /= 0 ) THEN
    308                    message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    309                             'res cloud_scheme = seifert_beheng'
    310                    CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    311303                ENDIF
    312304                unit = 'kg/kg'
     
    344336                            'res cloud_physics = .TRUE.'
    345337                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    346                 ELSEIF ( icloud_scheme /= 0 ) THEN
     338                ELSEIF ( .NOT. microphysics_seifert ) THEN
    347339                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    348340                            'res cloud_scheme = seifert_beheng'
    349341                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    350                 ELSEIF ( .NOT. precipitation )  THEN
    351                    message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    352                                     'res precipitation = .TRUE.'
    353                    CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    354342                ENDIF
    355343                unit = 'kg/kg'
  • palm/trunk/SOURCE/interaction_droplets_ptq.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Unused variables removed.
    2222!
    2323! Former revisions:
     
    127127
    128128       USE indices,                                                            &
    129            ONLY:  nxl, nxr, nyn, nys, nzb_2d, nzt
     129           ONLY:  nzb_2d, nzt
    130130
    131131       USE kinds,                                                              &
  • palm/trunk/SOURCE/lpm.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Tails removed.
     22!
     23! Initialization of sgs model not necessary for the use of cloud_droplets and
     24! use_sgs_for_particles.
     25!
     26! lpm_release_set integrated.
     27!
     28! Unused variabled removed.
    2229!
    2330! Former revisions:
     
    108115        ONLY:  lpm_exchange_horiz, lpm_move_particle
    109116
     117    USE lpm_init_mod,                                                          &
     118        ONLY: lpm_create_particle, PHASE_RELEASE
     119
    110120    USE lpm_pack_arrays_mod,                                                   &
    111121        ONLY:  lpm_pack_all_arrays
    112122
    113123    USE particle_attributes,                                                   &
    114         ONLY:  collision_kernel, deleted_particles, deleted_tails,             &
     124        ONLY:  collision_kernel, deleted_particles,                            &
    115125               dt_write_particle_data, dt_prel, end_time_prel,                 &
    116126               grid_particles, number_of_particles, number_of_particle_groups, &
    117127               particles, particle_groups, prt_count, trlp_count_sum,          &
    118                tail_mask, time_prel, time_sort_particles,                      &
     128               time_prel,                                                      &
    119129               time_write_particle_data, trlp_count_recv_sum, trnp_count_sum,  &
    120130               trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &
    121                trsp_count_sum, trsp_count_recv_sum, use_particle_tails,        &
     131               trsp_count_sum, trsp_count_recv_sum,                            &
    122132               use_sgs_for_particles, write_particle_statistics
    123133
     
    159169
    160170!
    161 !-- Initialize arrays for marking those particles/tails to be deleted after the
     171!-- Initialize arrays for marking those particles to be deleted after the
    162172!-- (sub-) timestep
    163173    deleted_particles = 0
    164 
    165     IF ( use_particle_tails )  THEN
    166        tail_mask = .TRUE.
    167     ENDIF
    168     deleted_tails = 0
    169 
    170174
    171175!
     
    202206    IF ( time_prel >= dt_prel  .AND.  end_time_prel > simulated_time )  THEN
    203207
    204        CALL lpm_release_set
     208       CALL lpm_create_particle(PHASE_RELEASE)
    205209!
    206210!--    The MOD function allows for changes in the output interval with
     
    240244!--    velocity variances)
    241245
    242        IF ( use_sgs_for_particles )  CALL lpm_init_sgs_tke
     246       IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
     247          CALL lpm_init_sgs_tke
     248       ENDIF
    243249
    244250       DO  i = nxl, nxr
     
    361367       deleted_particles = 0
    362368
    363        IF ( use_particle_tails )  THEN
    364           tail_mask     = .TRUE.
    365        ENDIF
    366        deleted_tails = 0
    367 
    368 
    369369       IF ( dt_3d_reached )  EXIT
    370370
     
    387387    IF ( collision_kernel == 'none' )  CALL lpm_set_attributes
    388388
    389 
    390389!
    391390!-- Set particle attributes defined by the user
    392391    CALL user_lpm_set_attributes
    393 
    394 
    395 !
    396 !-- particle tails currently not available
    397 !
    398 !-- If required, add the current particle positions to the particle tails
    399 !   IF ( use_particle_tails )  CALL lpm_extend_tails
    400 
    401392
    402393!
  • palm/trunk/SOURCE/lpm_advec.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Random velocity fluctuations for particles added. Terminal fall velocity
     22! for droplets is calculated from a parameterization (which is better than
     23! the previous, physically correct calculation, which demands a very short
     24! time step that is not used in the model).
     25!
     26! Unused variables deleted.
    2227!
    2328! Former revisions:
     
    7883
    7984    USE arrays_3d,                                                             &
    80         ONLY:  de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu,  &
    81                zw
     85        ONLY:  de_dx, de_dy, de_dz, diss, e, km, u, us, usws, v, vsws, w, zu, zw
    8286
    8387    USE cpulog
     
    8791    USE control_parameters,                                                    &
    8892        ONLY:  atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d,   &
    89                dt_3d_reached_l, dz, g, kappa, molecular_viscosity, topography, &
    90                u_gtrans, v_gtrans, simulated_time
     93               dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
    9194
    9295    USE grid_variables,                                                        &
     
    99102   
    100103    USE particle_attributes,                                                   &
    101         ONLY:  block_offset, c_0, density_ratio, dt_min_part, grid_particles,  &
     104        ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
    102105               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
    103106               particles, particle_groups, offset_ocean_nzt, sgs_wfu_part,     &
     
    140143    REAL(wp) ::  de_dz_int_l        !<
    141144    REAL(wp) ::  de_dz_int_u        !<
     145    REAL(wp) ::  diameter           !< diamter of droplet
    142146    REAL(wp) ::  diss_int_l         !<
    143147    REAL(wp) ::  diss_int_u         !<
     
    150154    REAL(wp) ::  exp_term           !<
    151155    REAL(wp) ::  gg                 !<
    152     REAL(wp) ::  height_int         !<
    153156    REAL(wp) ::  height_p           !<
    154     REAL(wp) ::  lagr_timescale     !<
     157    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
    155158    REAL(wp) ::  location(1:30,1:3) !<
    156159    REAL(wp) ::  random_gauss       !<
     160    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
     161    REAL(wp) ::  rg1                !< Gaussian distributed random number
     162    REAL(wp) ::  rg2                !< Gaussian distributed random number
     163    REAL(wp) ::  rg3                !< Gaussian distributed random number
     164    REAL(wp) ::  sigma              !< velocity standard deviation
    157165    REAL(wp) ::  u_int_l            !<
    158166    REAL(wp) ::  u_int_u            !<
     
    163171    REAL(wp) ::  w_int_l            !<
    164172    REAL(wp) ::  w_int_u            !<
     173    REAL(wp) ::  w_s                !< terminal velocity of droplets
    165174    REAL(wp) ::  x                  !<
    166175    REAL(wp) ::  y                  !<
    167     REAL(wp) ::  z_p                !<   
     176    REAL(wp) ::  z_p                !<
     177
     178    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
     179    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
     180    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
     181    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
     182    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
     183    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
    168184
    169185    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !<
     
    388404!-- Interpolate and calculate quantities needed for calculating the SGS
    389405!-- velocities
    390     IF ( use_sgs_for_particles )  THEN
     406    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    391407
    392408       IF ( topography == 'flat' )  THEN
     
    13361352!--          Update of the particle velocity
    13371353             IF ( cloud_droplets )  THEN
    1338                 exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /      &
    1339                            ( particles(n)%radius )**2 *                        &
    1340                            ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius &
    1341                              * SQRT( ( u_int(n) - particles(n)%speed_x )**2 +  &
    1342                                      ( v_int(n) - particles(n)%speed_y )**2 +  &
    1343                                      ( w_int(n) - particles(n)%speed_z )**2 )  &
    1344                                           / molecular_viscosity )**0.687_wp    &
    1345                            )
    1346 
    1347                 exp_term = EXP( -exp_arg * dt_particle(n) )
    1348              ELSEIF ( use_sgs_for_particles )  THEN
     1354!
     1355!--             Terminal velocity is computed for vertical direction (Rogers et
     1356!--             al., 1993, J. Appl. Meteorol.)
     1357                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
     1358                IF ( diameter <= d0_rog )  THEN
     1359                   w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
     1360                ELSE
     1361                   w_s = a_rog - b_rog * EXP( -c_rog * diameter )
     1362                ENDIF
     1363
     1364!
     1365!--             If selected, add random velocities following Soelch and Kaercher
     1366!--             (2010, Q. J. R. Meteorol. Soc.)
     1367                IF ( use_sgs_for_particles )  THEN
     1368                   lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
     1369                   RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
     1370                   sigma          = SQRT( e(kp,jp,ip) )
     1371
     1372                   rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1373                   rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1374                   rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1375
     1376                   particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
     1377                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg1
     1378                   particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
     1379                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg2
     1380                   particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
     1381                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg3
     1382
     1383                   particles(n)%speed_x = u_int(n) + particles(n)%rvar1
     1384                   particles(n)%speed_y = v_int(n) + particles(n)%rvar2
     1385                   particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
     1386                ELSE
     1387                   particles(n)%speed_x = u_int(n)
     1388                   particles(n)%speed_y = v_int(n)
     1389                   particles(n)%speed_z = w_int(n) - w_s
     1390                ENDIF
     1391
     1392             ELSE
     1393
     1394                IF ( use_sgs_for_particles )  THEN
     1395                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
     1396                   exp_term = EXP( -exp_arg * dt_particle(n) )
     1397                ELSE
     1398                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
     1399                   exp_term = particle_groups(particles(n)%group)%exp_term
     1400                ENDIF
     1401                particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
     1402                                       u_int(n) * ( 1.0_wp - exp_term )
     1403                particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
     1404                                       v_int(n) * ( 1.0_wp - exp_term )
     1405                particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
     1406                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
     1407                                       g / exp_arg ) * ( 1.0_wp - exp_term )
     1408             ENDIF
     1409
     1410          ENDIF
     1411
     1412       ENDDO
     1413   
     1414    ELSE
     1415
     1416       DO  n = 1, number_of_particles
     1417
     1418!--       Transport of particles with inertia
     1419          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
     1420          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
     1421          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
     1422!
     1423!--       Update of the particle velocity
     1424          IF ( cloud_droplets )  THEN
     1425!
     1426!--          Terminal velocity is computed for vertical direction (Rogers et al.,
     1427!--          1993, J. Appl. Meteorol.)
     1428             diameter = particles(n)%radius * 2000.0_wp !diameter in mm
     1429             IF ( diameter <= d0_rog )  THEN
     1430                w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
     1431             ELSE
     1432                w_s = a_rog - b_rog * EXP( -c_rog * diameter )
     1433             ENDIF
     1434
     1435!
     1436!--          If selected, add random velocities following Soelch and Kaercher
     1437!--          (2010, Q. J. R. Meteorol. Soc.)
     1438             IF ( use_sgs_for_particles )  THEN
     1439                 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
     1440                 RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
     1441                 sigma          = SQRT( e(kp,jp,ip) )
     1442
     1443                 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1444                 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1445                 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1446
     1447                 particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
     1448                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg1
     1449                 particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
     1450                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg2
     1451                 particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
     1452                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg3
     1453
     1454                 particles(n)%speed_x = u_int(n) + particles(n)%rvar1
     1455                 particles(n)%speed_y = v_int(n) + particles(n)%rvar2
     1456                 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
     1457             ELSE
     1458                 particles(n)%speed_x = u_int(n)
     1459                 particles(n)%speed_y = v_int(n)
     1460                 particles(n)%speed_z = w_int(n) - w_s
     1461             ENDIF
     1462
     1463          ELSE
     1464
     1465             IF ( use_sgs_for_particles )  THEN
    13491466                exp_arg  = particle_groups(particles(n)%group)%exp_arg
    13501467                exp_term = EXP( -exp_arg * dt_particle(n) )
     
    13531470                exp_term = particle_groups(particles(n)%group)%exp_term
    13541471             ENDIF
    1355              particles(n)%speed_x = particles(n)%speed_x * exp_term +          &
     1472             particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
    13561473                                    u_int(n) * ( 1.0_wp - exp_term )
    1357              particles(n)%speed_y = particles(n)%speed_y * exp_term +          &
     1474             particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
    13581475                                    v_int(n) * ( 1.0_wp - exp_term )
    1359              particles(n)%speed_z = particles(n)%speed_z * exp_term +          &
    1360                                     ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
    1361                                     g / exp_arg ) * ( 1.0_wp - exp_term )
     1476             particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
     1477                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
     1478                                    exp_arg ) * ( 1.0_wp - exp_term )
    13621479          ENDIF
    13631480
    1364        ENDDO
    1365    
    1366     ELSE
    1367 
    1368        DO  n = 1, number_of_particles
    1369 
    1370 !--       Transport of particles with inertia
    1371           particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
    1372           particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
    1373           particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
    1374 !
    1375 !--       Update of the particle velocity
    1376           IF ( cloud_droplets )  THEN
    1377 
    1378              exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /         &
    1379                         ( particles(n)%radius )**2 *                           &
    1380                         ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius *  &
    1381                           SQRT( ( u_int(n) - particles(n)%speed_x )**2 +       &
    1382                                 ( v_int(n) - particles(n)%speed_y )**2 +       &
    1383                                 ( w_int(n) - particles(n)%speed_z )**2 ) /     &
    1384                                          molecular_viscosity )**0.687_wp       &
    1385                         )
    1386 
    1387              exp_term = EXP( -exp_arg * dt_particle(n) )
    1388           ELSEIF ( use_sgs_for_particles )  THEN
    1389              exp_arg  = particle_groups(particles(n)%group)%exp_arg
    1390              exp_term = EXP( -exp_arg * dt_particle(n) )
    1391           ELSE
    1392              exp_arg  = particle_groups(particles(n)%group)%exp_arg
    1393              exp_term = particle_groups(particles(n)%group)%exp_term
    1394           ENDIF
    1395           particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
    1396                                  u_int(n) * ( 1.0_wp - exp_term )
    1397           particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
    1398                                  v_int(n) * ( 1.0_wp - exp_term )
    1399           particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
    1400                                  ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
    1401                                  exp_arg ) * ( 1.0_wp - exp_term )
    14021481       ENDDO
    14031482
  • palm/trunk/SOURCE/lpm_boundary_conds.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Tails removed. Unused variables removed.
    2222!
    2323! Former revisions:
     
    7575
    7676    USE control_parameters,                                                    &
    77         ONLY:  dz, message_string, particle_maximum_age, simulated_time
     77        ONLY:  dz, message_string, particle_maximum_age
    7878
    7979    USE cpulog,                                                                &
     
    8989
    9090    USE particle_attributes,                                                   &
    91         ONLY:  deleted_particles, deleted_tails, ibc_par_b, ibc_par_t,         &
    92                number_of_particles, particles,                                 &
    93                particle_tail_coordinates, particle_type, offset_ocean_nzt_m1,  &
    94                tail_mask, use_particle_tails, use_sgs_for_particles
     91        ONLY:  deleted_particles, ibc_par_b, ibc_par_t, number_of_particles,   &
     92               particles, particle_type, offset_ocean_nzt_m1,                  &
     93               use_sgs_for_particles
    9594
    9695    USE pegrid
     
    119118    INTEGER(iwp) ::  k5             !<
    120119    INTEGER(iwp) ::  n              !<
    121     INTEGER(iwp) ::  nn             !<
    122120    INTEGER(iwp) ::  t_index        !<
    123121    INTEGER(iwp) ::  t_index_number !<
     
    150148       DO  n = 1, number_of_particles
    151149
    152           nn = particles(n)%tail_id
    153 
    154150!
    155151!--       Stop if particles have moved further than the length of one
     
    171167             particles(n)%particle_mask  = .FALSE.
    172168             deleted_particles = deleted_particles + 1
    173              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    174                 tail_mask(nn) = .FALSE.
    175                 deleted_tails = deleted_tails + 1
    176              ENDIF
    177169          ENDIF
    178170
     
    183175                particles(n)%particle_mask  = .FALSE.
    184176                deleted_particles = deleted_particles + 1
    185                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    186                    tail_mask(nn) = .FALSE.
    187                    deleted_tails = deleted_tails + 1
    188                 ENDIF
    189177             ELSEIF ( ibc_par_t == 2 )  THEN
    190178!
     
    196184                   particles(n)%rvar3 = -particles(n)%rvar3
    197185                ENDIF
    198                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    199                    particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
    200                                                particle_tail_coordinates(1,3,nn)
    201                 ENDIF
    202186             ENDIF
    203187          ENDIF
     
    209193                particles(n)%particle_mask  = .FALSE.
    210194                deleted_particles = deleted_particles + 1
    211                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    212                    tail_mask(nn) = .FALSE.
    213                    deleted_tails = deleted_tails + 1
    214                 ENDIF
    215195             ELSEIF ( ibc_par_b == 2 )  THEN
    216196!
     
    221201                     particles(n)%rvar3 < 0.0_wp )  THEN
    222202                   particles(n)%rvar3 = -particles(n)%rvar3
    223                 ENDIF
    224                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    225                    particle_tail_coordinates(1,3,nn) = 2.0_wp * zu(nz) - &
    226                                                particle_tail_coordinates(1,3,nn)
    227                 ENDIF
    228                 IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    229                    particle_tail_coordinates(1,3,nn) = 2.0_wp * zw(0) - &
    230                                                particle_tail_coordinates(1,3,nn)
    231203                ENDIF
    232204             ENDIF
  • palm/trunk/SOURCE/lpm_calc_liquid_water_content.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Unused variables removed.
    2222!
    2323! Former revisions:
     
    8585    INTEGER(iwp) ::  k   !<
    8686    INTEGER(iwp) ::  n   !<
    87     INTEGER(iwp) ::  psi !<
    8887
    8988    CALL cpu_log( log_point_s(45), 'lpm_calc_ql', 'start' )
  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! PALM kernel has been deleted.
     22! Bugfix in the calculation of the turbulent enhancement factor of the
     23! collection efficiency.
     24!
     25! Unused variables removed.
    2226!
    2327! Former revisions:
     
    96100!> This module calculates collision efficiencies either due to pure gravitational
    97101!> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
    98 !> including the effects of (SGS) turbulence (Wang kernel, see Wang and
    99 !> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
     102!> including the effects of turbulence (Wang kernel, see Wang and
     103!> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8, and Ayala et al., 2008:
     104!> New J. Phys., 10, 075016). The original code has been
    100105!> provided by L.-P. Wang but is substantially reformatted and speed optimized
    101106!> here.
    102 !>
    103 !> @attention Physical quantities (like g, densities, etc.) used in this module
    104 !>            still have to be adjusted to those values used in the main PALM code.
    105 !>            Also, quantities in CGS-units should be converted to SI-units eventually.
    106107!------------------------------------------------------------------------------!
    107108 MODULE lpm_collision_kernels_mod
     
    114115
    115116    USE particle_attributes,                                                   &
    116         ONLY:  collision_kernel, dissipation_classes, particles, radius_classes
     117        ONLY:  collision_kernel, dissipation_classes, particles,               &
     118               radius_classes
    117119
    118120    USE pegrid
     
    123125    PRIVATE
    124126
    125     PUBLIC  ckernel, collision_efficiency_rogers, init_kernels, &
    126             rclass_lbound, rclass_ubound, recalculate_kernel
     127    PUBLIC  ckernel, init_kernels, rclass_lbound, rclass_ubound,              &
     128            recalculate_kernel
    127129
    128130    REAL(wp) ::  epsilon       !<
    129     REAL(wp) ::  eps2          !<
    130131    REAL(wp) ::  rclass_lbound !<
    131132    REAL(wp) ::  rclass_ubound !<
    132133    REAL(wp) ::  urms          !<
    133     REAL(wp) ::  urms2         !<
    134 
    135     REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass !<
    136     REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass !<
    137     REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf     !<
     134
     135    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass  !< dissipation rate class
     136    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass  !< radius class
     137    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf      !<
    138138   
    139     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ec       !<
    140     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ecf      !<
    141     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gck      !<
    142     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hkernel  !<
    143     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hwratio  !<
     139    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ec        !<
     140    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ecf       !<
     141    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gck       !<
     142    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hkernel   !<
     143    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hwratio   !<
    144144   
    145     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  ckernel !<
     145    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
    146146
    147147    SAVE
     
    149149!
    150150!-- Public interfaces
    151     INTERFACE collision_efficiency_rogers
    152        MODULE PROCEDURE collision_efficiency_rogers
    153     END INTERFACE collision_efficiency_rogers
    154 
    155151    INTERFACE init_kernels
    156152       MODULE PROCEDURE init_kernels
     
    186182       IF ( collision_kernel(6:9) == 'fast' )  THEN
    187183
    188           ALLOCATE( ckernel(1:radius_classes,1:radius_classes,               &
    189                     0:dissipation_classes), epsclass(1:dissipation_classes), &
     184          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,                 &
     185                    0:dissipation_classes), epsclass(1:dissipation_classes),   &
    190186                    radclass(1:radius_classes) )
    191187
     
    195191          rclass_lbound = LOG( 1.0E-6_wp )
    196192          rclass_ubound = LOG( 2.0E-4_wp )
    197           radclass(1)   = 1.0E-6_wp
     193          radclass(1)   = EXP( rclass_lbound )
    198194          DO  i = 2, radius_classes
    199195             radclass(i) = EXP( rclass_lbound +                                &
     
    377373! Description:
    378374! ------------
    379 !> Calculation of gck
    380 !> This is from Aayala 2008b, page 37ff.
    381 !> Necessary input parameters: water density, radii of droplets, air density,
    382 !> air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
    383 !> gravitational acceleration  --> to be replaced by PALM parameters
     375!> Calculation of effects of turbulence on the geometric collision kernel
     376!> (by including the droplets' average radial relative velocities and their
     377!> radial distribution function) following the analytic model by Aayala et al.
     378!> (2008, New J. Phys.). For details check the second part 2 of the publication,
     379!> page 37ff.
     380!>
     381!> Input parameters, which need to be replaced by PALM parameters:
     382!>    water density, air density
    384383!------------------------------------------------------------------------------!
    385384    SUBROUTINE turbsd
     
    392391
    393392       IMPLICIT NONE
    394        
    395        LOGICAL, SAVE ::  first = .TRUE. !<
    396393
    397394       INTEGER(iwp) ::  i     !<
     
    442439       REAL(wp) ::  z         !<
    443440
    444        REAL(wp), DIMENSION(1:radius_classes) ::  st  !<
    445        REAL(wp), DIMENSION(1:radius_classes) ::  tau !<
    446        
    447 !
    448 !--    Initial assignment of constants
    449        IF ( first )  THEN
    450 
    451           first = .FALSE.
    452 
    453        ENDIF
    454 
    455        lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon ) ! in m
     441       REAL(wp), DIMENSION(1:radius_classes) ::  st  !< Stokes number
     442       REAL(wp), DIMENSION(1:radius_classes) ::  tau !< inertial time scale
     443
     444       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon )
    456445       lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
    457        tl        = urms**2 / epsilon                       ! in s
    458        lf        = 0.5_wp * urms**3 / epsilon              ! in m
    459        tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
    460        eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp          ! in m
     446       tl        = urms**2 / epsilon
     447       lf        = 0.5_wp * urms**3 / epsilon
     448       tauk      = SQRT( molecular_viscosity / epsilon )
     449       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp
    461450       vk        = eta / tauk
    462451
    463452       ao = ( 11.0_wp + 7.0_wp * lambda_re ) / ( 205.0_wp + lambda_re )
    464        tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk  ! in s
    465 
    466        CALL fallg    ! gives winf in m/s
     453       tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk
     454
     455!
     456!--    Get terminal velocity of droplets
     457       CALL fallg
    467458
    468459       DO  i = 1, radius_classes
    469           tau(i) = winf(i) / g    ! in s
    470           st(i)  = tau(i) / tauk
     460          tau(i) = winf(i) / g    ! inertial time scale
     461          st(i)  = tau(i) / tauk  ! Stokes number
    471462       ENDDO
    472463
    473464!
    474 !--    Calculate wr (from Aayala 2008b, page 38f)
     465!--    Calculate average radial relative velocity at contact (wrfin)
    475466       z   = tt / tl
    476467       be  = SQRT( 2.0_wp ) * lambda / lf
    477468       bbb = SQRT( 1.0_wp - 2.0_wp * be**2 )
    478469       d1  = ( 1.0_wp + bbb ) / ( 2.0_wp * bbb )
    479        e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp   ! in m
     470       e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp
    480471       d2  = ( 1.0_wp - bbb ) * 0.5_wp / bbb
    481        e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp   ! in m
     472       e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp
    482473       ccc = SQRT( 1.0_wp - 2.0_wp * z**2 )
    483474       b1  = ( 1.0_wp + ccc ) * 0.5_wp / ccc
    484        c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp   ! in s
     475       c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp
    485476       b2  = ( 1.0_wp - ccc ) * 0.5_wp / ccc
    486        c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp   ! in s
     477       c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp
    487478
    488479       DO  i = 1, radius_classes
    489480
    490           v1 = winf(i)        ! in m/s
    491           t1 = tau(i)         ! in s
     481          v1 = winf(i)
     482          t1 = tau(i)
    492483
    493484          DO  j = 1, i
    494485             rrp = radclass(i) + radclass(j)
    495              v2  = winf(j)                                 ! in m/s
    496              t2  = tau(j)                                  ! in s
     486             v2  = winf(j)
     487             t2  = tau(j)
    497488
    498489             v1xysq  = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) &
    499490                     - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1)
    500              v1xysq  = v1xysq * urms**2 / t1                ! in m**2/s**2
    501              vrms1xy = SQRT( v1xysq )                       ! in m/s
     491             v1xysq  = v1xysq * urms**2 / t1
     492             vrms1xy = SQRT( v1xysq )
    502493
    503494             v2xysq  = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) &
    504495                     - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2)
    505              v2xysq  = v2xysq * urms**2 / t2                ! in m**2/s**2
    506              vrms2xy = SQRT( v2xysq )                       ! in m/s
     496             v2xysq  = v2xysq * urms**2 / t2
     497             vrms2xy = SQRT( v2xysq )
    507498
    508499             IF ( winf(i) >= winf(j) )  THEN
     
    523514                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
    524515             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
    525              v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)   ! in m**2/s**2
    526              wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy ! in m**2/s**2
     516             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)
     517             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy
    527518             IF ( wrtur2xy < 0.0_wp )  wrtur2xy = 0.0_wp
    528519             wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
    529              wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) ) ! in m/s
    530 
    531 !
    532 !--          Calculate gr
     520             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) )
     521
     522!
     523!--          Calculate radial distribution function (grfin)
    533524             IF ( st(j) > st(i) )  THEN
    534525                sst = st(j)
     
    546537             ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
    547538             fao_gr = 20.115_wp * SQRT( ao_gr / lambda_re )
    548              rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
     539             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta
    549540
    550541             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5_wp )
    551542             IF ( grfin < 1.0_wp )  grfin = 1.0_wp
    552543
    553              gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin        ! in cm**3/s
     544!
     545!--          Calculate general collection kernel (without the consideration of
     546!--          collection efficiencies)
     547             gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin
    554548             gck(j,i) = gck(i,j)
    555549
     
    559553    END SUBROUTINE turbsd
    560554
    561 
    562 !------------------------------------------------------------------------------!
    563 ! Description:
    564 ! ------------
    565 !> phi_w as a function
    566 !------------------------------------------------------------------------------!
    567555    REAL(wp) FUNCTION phi_w( a, b, vsett, tau0 )
    568 
     556!
     557!--    Function used in the Ayala et al. (2008) analytical model for turbulent
     558!--    effects on the collision kernel
    569559       IMPLICIT NONE
    570560
     
    576566
    577567       aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b
    578        phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2  ! in s
     568       phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2
    579569
    580570    END FUNCTION phi_w
    581571
    582 
    583 !------------------------------------------------------------------------------!
    584 ! Description:
    585 ! ------------
    586 !> zhi as a function
    587 !------------------------------------------------------------------------------!
    588572    REAL(wp) FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
    589 
     573!
     574!--    Function used in the Ayala et al. (2008) analytical model for turbulent
     575!--    effects on the collision kernel
    590576       IMPLICIT NONE
    591577
     
    613599             b / aa3**2 + ( 4.0_wp / aa4 - 1.0_wp / aa5**2 - 1.0_wp / aa1**2 ) &
    614600             * vsett2 * 0.5_wp / b /aa6 + ( 2.0_wp * ( b / aa2 - b / aa1 ) -   &
    615              vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3    ! in s**2
     601             vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3
    616602
    617603    END FUNCTION zhi
     
    621607! Description:
    622608! ------------
    623 !> Calculation of terminal velocity winf following Equations 10-138 to 10-145
    624 !> from (Pruppacher and Klett, 1997)
     609!> Parameterization of terminal velocity following Rogers et al. (1993, J. Appl.
     610!> Meteorol.)
    625611!------------------------------------------------------------------------------!
    626612    SUBROUTINE fallg
    627  
    628        USE cloud_parameters,                                                   &
    629            ONLY:  rho_l
    630    
    631        USE control_parameters,                                                 &
    632            ONLY:  g
    633613
    634614       USE particle_attributes,                                                &
    635615           ONLY:  radius_classes
    636616
    637 
    638617       IMPLICIT NONE
    639618
    640        INTEGER(iwp) ::  i !<
    641        INTEGER(iwp) ::  j !<
    642 
    643        LOGICAL, SAVE ::  first = .TRUE. !<
    644 
    645        REAL(wp), SAVE ::  cunh  !<
    646        REAL(wp), SAVE ::  eta   !<
    647        REAL(wp), SAVE ::  phy   !<
    648        REAL(wp), SAVE ::  py    !<
    649        REAL(wp), SAVE ::  rho_a !<
    650        REAL(wp), SAVE ::  sigma !<
    651        REAL(wp), SAVE ::  stb   !<
    652        REAL(wp), SAVE ::  stok  !<
    653        REAL(wp), SAVE ::  xlamb !<
    654 
    655        REAL(wp) ::  bond        !<
    656        REAL(wp) ::  x           !<
    657        REAL(wp) ::  xrey        !<
    658        REAL(wp) ::  y           !<
    659 
    660        REAL(wp), DIMENSION(1:7), SAVE  ::  b !<
    661        REAL(wp), DIMENSION(1:6), SAVE  ::  c !<
    662 
    663 !
    664 !--    Initial assignment of constants
    665        IF ( first )  THEN
    666 
    667           first = .FALSE.
    668           b = (/  -0.318657E1_wp,   0.992696E0_wp,  -0.153193E-2_wp, &
    669                   -0.987059E-3_wp, -0.578878E-3_wp,  0.855176E-4_wp, &
    670                   -0.327815E-5_wp /)
    671           c = (/  -0.500015E1_wp,   0.523778E1_wp,  -0.204914E1_wp,   &
    672                    0.475294E0_wp,  -0.542819E-1_wp,  0.238449E-2_wp /)
    673 
    674 !
    675 !--       Parameter values for p = 1013,25 hPa and T = 293,15 K
    676           eta   = 1.818E-5_wp         ! in kg/(m s)
    677           xlamb = 6.6E-8_wp           ! in m
    678           rho_a = 1.204_wp            ! in kg/m**3
    679           cunh  = 1.26_wp * xlamb     ! in m
    680           sigma = 0.07363_wp          ! in kg/s**2
    681           stok  = 2.0_wp  * g * ( rho_l - rho_a ) / ( 9.0_wp * eta ) ! in 1/(m s)
    682           stb   = 32.0_wp * rho_a * ( rho_l - rho_a) * g / (3.0_wp * eta * eta)
    683           phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
    684           py    = phy**( 1.0_wp / 6.0_wp )
    685 
    686        ENDIF
     619       INTEGER(iwp) ::  j                            !<
     620
     621       REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp    !< parameter
     622       REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp   !< parameter
     623       REAL(wp), PARAMETER ::  a_rog     = 9.65_wp   !< parameter
     624       REAL(wp), PARAMETER ::  b_rog     = 10.43_wp  !< parameter
     625       REAL(wp), PARAMETER ::  c_rog     = 0.6_wp    !< parameter
     626       REAL(wp), PARAMETER ::  d0_rog    = 0.745_wp  !< seperation diameter
     627
     628       REAL(wp)            ::  diameter              !< droplet diameter in mm
     629
    687630
    688631       DO  j = 1, radius_classes
    689632
    690           IF ( radclass(j) <= 1.0E-5_wp )  THEN
    691 
    692              winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) )
    693 
    694           ELSEIF ( radclass(j) > 1.0E-5_wp  .AND.  radclass(j) <= 5.35E-4_wp )  THEN
    695 
    696              x = LOG( stb * radclass(j)**3 )
    697              y = 0.0_wp
    698 
    699              DO  i = 1, 7
    700                 y = y + b(i) * x**(i-1)
    701              ENDDO
    702 !
    703 !--          Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418)
    704 !--          for correct version see (Beard, 1976)
    705              xrey = ( 1.0_wp + cunh / radclass(j) ) * EXP( y )
    706 
    707              winf(j) = xrey * eta / ( 2.0_wp * rho_a * radclass(j) )
    708 
    709           ELSEIF ( radclass(j) > 5.35E-4_wp )  THEN
    710 
    711              IF ( radclass(j) > 0.0035_wp )  THEN
    712                 bond = g * ( rho_l - rho_a ) * 0.0035_wp**2 / sigma
    713              ELSE
    714                bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma
    715              ENDIF
    716 
    717              x = LOG( 16.0_wp * bond * py / 3.0_wp )
    718              y = 0.0_wp
    719 
    720              DO  i = 1, 6
    721                 y = y + c(i) * x**(i-1)
    722              ENDDO
    723 
    724              xrey = py * EXP( y )
    725 
    726              IF ( radclass(j) > 0.0035_wp )  THEN
    727                 winf(j) = xrey * eta / ( 2.0_wp * rho_a * 0.0035_wp )
    728              ELSE
    729                 winf(j) = xrey * eta / ( 2.0_wp * rho_a * radclass(j) )
    730              ENDIF
    731 
     633          diameter = radclass(j) * 2000.0_wp
     634
     635          IF ( diameter <= d0_rog )  THEN
     636             winf(j) = k_cap_rog * diameter * ( 1.0_wp -                       &
     637                                                EXP( -k_low_rog * diameter ) )
     638          ELSE
     639             winf(j) = a_rog - b_rog * EXP( -c_rog * diameter )
    732640          ENDIF
    733641
     
    740648! Description:
    741649! ------------
    742 !> Calculation of collision efficiencies for the Hall kernel
     650!> Interpolation of collision efficiencies (Hall, 1980, J. Atmos. Sci.)
    743651!------------------------------------------------------------------------------!
    744652    SUBROUTINE effic
     
    776684
    777685         first = .FALSE.
    778          r0  = (/   6.0_wp,   8.0_wp,  10.0_wp, 15.0_wp,  20.0_wp,  25.0_wp,  &
    779                    30.0_wp,  40.0_wp,  50.0_wp, 60.0_wp,  70.0_wp, 100.0_wp,  &
     686         r0  = (/   6.0_wp,   8.0_wp,  10.0_wp, 15.0_wp,  20.0_wp,  25.0_wp,   &
     687                   30.0_wp,  40.0_wp,  50.0_wp, 60.0_wp,  70.0_wp, 100.0_wp,   &
    780688                  150.0_wp, 200.0_wp, 300.0_wp /)
    781689
    782          rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp,       &
    783                   0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp,       &
    784                   0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp,       &
     690         rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp,        &
     691                  0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp,        &
     692                  0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp,        &
    785693                  0.90_wp, 0.95_wp, 1.00_wp /)
    786694
    787          ecoll(:,1)  = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, &
    788                           0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, &
     695         ecoll(:,1)  = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp,    &
     696                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp,    &
    789697                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp /)
    790          ecoll(:,2)  = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp, &
    791                           0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp, &
     698         ecoll(:,2)  = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp,    &
     699                          0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp,    &
    792700                          0.200_wp, 0.500_wp, 0.770_wp, 0.870_wp, 0.970_wp /)
    793          ecoll(:,3)  = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp, &
    794                           0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp, &
     701         ecoll(:,3)  = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp,    &
     702                          0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp,    &
    795703                          0.580_wp, 0.790_wp, 0.930_wp, 0.960_wp, 1.000_wp /)
    796          ecoll(:,4)  = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp, &
    797                           0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp, &
     704         ecoll(:,4)  = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp,    &
     705                          0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp,    &
    798706                          0.750_wp, 0.910_wp, 0.970_wp, 0.980_wp, 1.000_wp /)
    799          ecoll(:,5)  = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp, &
    800                           0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp, &
     707         ecoll(:,5)  = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp,    &
     708                          0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp,    &
    801709                          0.840_wp, 0.950_wp, 0.970_wp, 1.000_wp, 1.000_wp /)
    802          ecoll(:,6)  = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp, &
    803                           0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp, &
     710         ecoll(:,6)  = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp,    &
     711                          0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp,    &
    804712                          0.880_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    805          ecoll(:,7)  = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp, &
    806                           0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp, &
     713         ecoll(:,7)  = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp,    &
     714                          0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp,    &
    807715                          0.900_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    808          ecoll(:,8)  = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp, &
    809                           0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp, &
     716         ecoll(:,8)  = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp,    &
     717                          0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp,    &
    810718                          0.920_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    811          ecoll(:,9)  = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp, &
    812                           0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp, &
     719         ecoll(:,9)  = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp,    &
     720                          0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp,    &
    813721                          0.940_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    814          ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp, &
    815                           0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
     722         ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp,    &
     723                          0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
    816724                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    817          ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp, &
    818                           0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
     725         ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp,    &
     726                          0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
    819727                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    820          ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp, &
    821                           0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp, &
     728         ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp,    &
     729                          0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
    822730                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    823          ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp, &
    824                           0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp, &
     731         ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp,    &
     732                          0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp,    &
    825733                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    826          ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp, &
    827                           0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp, &
     734         ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp,    &
     735                          0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp,    &
    828736                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    829          ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp, &
    830                           0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp, &
     737         ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp,    &
     738                          0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp,    &
    831739                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    832          ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp, &
    833                           0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp, &
     740         ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp,    &
     741                          0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp,    &
    834742                          0.970_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    835          ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp, &
    836                           0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp, &
     743         ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp,    &
     744                          0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp,    &
    837745                          1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
    838          ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp, &
    839                           0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp, &
     746         ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp,    &
     747                          0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp,    &
    840748                          1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp /)
    841          ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp, &
    842                           0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp, &
     749         ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp,    &
     750                          0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp,    &
    843751                          1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp /)
    844          ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, &
    845                           0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp, &
     752         ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp,    &
     753                          0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp,    &
    846754                          2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp /)
    847          ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, &
    848                           0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp, &
     755         ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp,    &
     756                          0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp,    &
    849757                          4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp /)
    850758       ENDIF
     
    852760!
    853761!--    Calculate the radius class index of particles with respect to array r
    854 !--    Radius has to be in µm
     762!--    Radius has to be in microns
    855763       ALLOCATE( ira(1:radius_classes) )
    856764       DO  j = 1, radius_classes
     
    867775!
    868776!--    Two-dimensional linear interpolation of the collision efficiency.
    869 !--    Radius has to be in µm
     777!--    Radius has to be in microns
    870778       DO  j = 1, radius_classes
    871779          DO  i = 1, j
     
    878786             IF ( ir < 16 )  THEN
    879787                IF ( ir >= 2 )  THEN
    880                    pp = ( ( MAX(radclass(j),radclass(i)) * 1.0E06_wp ) - r0(ir-1) ) / &
    881                         ( r0(ir) - r0(ir-1) )
     788                   pp = ( ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp ) -    &
     789                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
    882790                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    883791                   ec(j,i) = ( 1.0_wp - pp ) * ( 1.0_wp - qq )                 &
     
    911819! Description:
    912820! ------------
    913 !> Calculation of enhancement factor for collision efficencies due to turbulence
     821!> Interpolation of turbulent enhancement factor for collision efficencies
     822!> following Wang and Grabowski (2009, Atmos. Sci. Let.)
    914823!------------------------------------------------------------------------------!
    915824    SUBROUTINE turb_enhance_eff
     
    940849
    941850       REAL(wp), DIMENSION(1:11), SAVE ::  rat           !<
    942        
    943851       REAL(wp), DIMENSION(1:7), SAVE  ::  r0            !<
    944852       
     
    958866                   0.7_wp, 0.8_wp, 0.9_wp, 1.0_wp /)
    959867!
    960 !--       for 100 cm**2/s**3
     868!--       Tabulated turbulent enhancement factor at 100 cm**2/s**3
    961869          ecoll_100(:,1)  = (/  1.74_wp,   1.74_wp,   1.773_wp, 1.49_wp,  &
    962870                                1.207_wp,  1.207_wp,  1.0_wp /)
     
    979887          ecoll_100(:,10) = (/  1.570_wp,  1.570_wp,  1.244_wp, 1.166_wp, &
    980888                                1.088_wp,  1.088_wp,  1.0_wp /)
    981           ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,  8.61_wp,  &
     889          ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,   8.61_wp,  &
    982890                                2.60_wp,   2.60_wp,   1.0_wp /)
    983891!
    984 !--       for 400 cm**2/s**3
     892!--       Tabulated turbulent enhancement factor at 400 cm**2/s**3
    985893          ecoll_400(:,1)  = (/  4.976_wp,  4.976_wp,  3.593_wp,  2.519_wp, &
    986894                                1.445_wp,  1.445_wp,  1.0_wp /)
     
    1010918!
    1011919!--    Calculate the radius class index of particles with respect to array r0
    1012 !--    Radius has to be in µm
     920!--    The droplet radius has to be given in microns.
    1013921       ALLOCATE( ira(1:radius_classes) )
    1014922
     
    1025933
    1026934!
    1027 !--    Two-dimensional linear interpolation of the collision efficiencies
    1028 !--    Radius has to be in µm
     935!--    Two-dimensional linear interpolation of the turbulent enhancement factor.
     936!--    The droplet radius has to be given in microns.
    1029937       DO  j =  1, radius_classes
    1030938          DO  i = 1, j
     
    1040948             ENDDO
    1041949
    1042              y1 = 0.0001_wp      ! for 0 m**2/s**3
     950             y1 = 1.0_wp  ! turbulent enhancement factor at 0 m**2/s**3
    1043951
    1044952             IF ( ir < 8 )  THEN
    1045953                IF ( ir >= 2 )  THEN
    1046                    pp = ( MAX(radclass(j),radclass(i))*1.0E6_wp - r0(ir-1) ) / &
    1047                         ( r0(ir) - r0(ir-1) )
     954                   pp = ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp - &
     955                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
    1048956                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    1049957                   y2 = ( 1.0_wp - pp ) * ( 1.0_wp - qq ) * ecoll_100(ir-1,iq-1) + &
     
    1066974             ENDIF
    1067975!
    1068 !--          Linear interpolation of dissipation rate in m**2/s**3
     976!--          Linear interpolation of turbulent enhancement factor
    1069977             IF ( epsilon <= 0.01_wp )  THEN
    1070978                ecf(j,i) = ( epsilon - 0.01_wp ) / ( 0.0_wp  - 0.01_wp ) * y1 &
     
    1087995    END SUBROUTINE turb_enhance_eff
    1088996
    1089 
    1090 
    1091 !------------------------------------------------------------------------------!
    1092 ! Description:
    1093 ! ------------
    1094 !> Collision efficiencies from table 8.2 in Rogers and Yau (1989, 3rd edition).
    1095 !> Values are calculated from table by bilinear interpolation.
    1096 !------------------------------------------------------------------------------!
    1097  
    1098     SUBROUTINE collision_efficiency_rogers( mean_r, r, e)
    1099 
    1100        IMPLICIT NONE
    1101 
    1102        INTEGER(iwp)  ::  i !<
    1103        INTEGER(iwp)  ::  j !<
    1104        INTEGER(iwp)  ::  k !<
    1105 
    1106        LOGICAL, SAVE ::  first = .TRUE. !<
    1107 
    1108        REAL(wp)      ::  aa      !<
    1109        REAL(wp)      ::  bb      !<
    1110        REAL(wp)      ::  cc      !<
    1111        REAL(wp)      ::  dd      !<
    1112        REAL(wp)      ::  dx      !<
    1113        REAL(wp)      ::  dy      !<
    1114        REAL(wp)      ::  e       !<
    1115        REAL(wp)      ::  gg      !<
    1116        REAL(wp)      ::  mean_r  !<
    1117        REAL(wp)      ::  mean_rm !<
    1118        REAL(wp)      ::  r       !<
    1119        REAL(wp)      ::  rm      !<
    1120        REAL(wp)      ::  x       !<
    1121        REAL(wp)      ::  y       !<
    1122  
    1123        REAL(wp), DIMENSION(1:9), SAVE      ::  collected_r = 0.0_wp !<
    1124        
    1125        REAL(wp), DIMENSION(1:19), SAVE     ::  collector_r = 0.0_wp !<
    1126        
    1127        REAL(wp), DIMENSION(1:9,1:19), SAVE ::  ef = 0.0_wp          !<
    1128 
    1129        mean_rm = mean_r * 1.0E06_wp
    1130        rm      = r      * 1.0E06_wp
    1131 
    1132        IF ( first )  THEN
    1133 
    1134           collected_r = (/    2.0_wp,    3.0_wp,    4.0_wp,    6.0_wp,    8.0_wp, &
    1135                              10.0_wp,   15.0_wp,   20.0_wp,   25.0_wp /)
    1136           collector_r = (/   10.0_wp,   20.0_wp,   30.0_wp,   40.0_wp,   50.0_wp, &
    1137                              60.0_wp,   80.0_wp,  100.0_wp,  150.0_wp,  200.0_wp, &
    1138                             300.0_wp,  400.0_wp,  500.0_wp,  600.0_wp, 1000.0_wp, &
    1139                            1400.0_wp, 1800.0_wp, 2400.0_wp, 3000.0_wp /)
    1140 
    1141           ef(:,1)  = (/ 0.017_wp, 0.027_wp, 0.037_wp, 0.052_wp, 0.052_wp,      &
    1142                         0.052_wp, 0.052_wp, 0.0_wp,   0.0_wp /)
    1143           ef(:,2)  = (/ 0.001_wp, 0.016_wp, 0.027_wp, 0.060_wp, 0.12_wp,       &
    1144                         0.17_wp,  0.17_wp,  0.17_wp,  0.0_wp /)
    1145           ef(:,3)  = (/ 0.001_wp, 0.001_wp, 0.02_wp,  0.13_wp,  0.28_wp,       &
    1146                         0.37_wp,  0.54_wp,  0.55_wp,  0.47_wp/)
    1147           ef(:,4)  = (/ 0.001_wp, 0.001_wp, 0.02_wp,  0.23_wp,  0.4_wp,        &
    1148                         0.55_wp,  0.7_wp,   0.75_wp,  0.75_wp/)
    1149           ef(:,5)  = (/ 0.01_wp,  0.01_wp,  0.03_wp,  0.3_wp,   0.4_wp,        &
    1150                         0.58_wp,  0.73_wp,  0.75_wp,  0.79_wp/)
    1151           ef(:,6)  = (/ 0.01_wp,  0.01_wp,  0.13_wp,  0.38_wp,  0.57_wp,       &
    1152                         0.68_wp,  0.80_wp,  0.86_wp,  0.91_wp/)
    1153           ef(:,7)  = (/ 0.01_wp,  0.085_wp, 0.23_wp,  0.52_wp,  0.68_wp,       &
    1154                         0.76_wp,  0.86_wp,  0.92_wp,  0.95_wp/)
    1155           ef(:,8)  = (/ 0.01_wp,  0.14_wp,  0.32_wp,  0.60_wp,  0.73_wp,       &
    1156                         0.81_wp,  0.90_wp,  0.94_wp,  0.96_wp/)
    1157           ef(:,9)  = (/ 0.025_wp, 0.25_wp,  0.43_wp,  0.66_wp,  0.78_wp,       &
    1158                         0.83_wp,  0.92_wp,  0.95_wp,  0.96_wp/)
    1159           ef(:,10) = (/ 0.039_wp, 0.3_wp,   0.46_wp,  0.69_wp,  0.81_wp,       &
    1160                         0.87_wp,  0.93_wp,  0.95_wp,  0.96_wp/)
    1161           ef(:,11) = (/ 0.095_wp, 0.33_wp,  0.51_wp,  0.72_wp,  0.82_wp,       &
    1162                         0.87_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
    1163           ef(:,12) = (/ 0.098_wp, 0.36_wp,  0.51_wp,  0.73_wp,  0.83_wp,       &
    1164                         0.88_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
    1165           ef(:,13) = (/ 0.1_wp,   0.36_wp,  0.52_wp,  0.74_wp,  0.83_wp,       &
    1166                         0.88_wp,  0.93_wp,  0.96_wp,  0.97_wp/)
    1167           ef(:,14) = (/ 0.17_wp,  0.4_wp,   0.54_wp,  0.72_wp,  0.83_wp,       &
    1168                         0.88_wp,  0.94_wp,  0.98_wp,  1.0_wp /)
    1169           ef(:,15) = (/ 0.15_wp,  0.37_wp,  0.52_wp,  0.74_wp,  0.82_wp,       &
    1170                         0.88_wp,  0.94_wp,  0.98_wp,  1.0_wp /)
    1171           ef(:,16) = (/ 0.11_wp,  0.34_wp,  0.49_wp,  0.71_wp,  0.83_wp,       &
    1172                         0.88_wp,  0.94_wp,  0.95_wp,  1.0_wp /)
    1173           ef(:,17) = (/ 0.08_wp,  0.29_wp,  0.45_wp,  0.68_wp,  0.8_wp,        &
    1174                         0.86_wp,  0.96_wp,  0.94_wp,  1.0_wp /)
    1175           ef(:,18) = (/ 0.04_wp,  0.22_wp,  0.39_wp,  0.62_wp,  0.75_wp,       &
    1176                         0.83_wp,  0.92_wp,  0.96_wp,  1.0_wp /)
    1177           ef(:,19) = (/ 0.02_wp,  0.16_wp,  0.33_wp,  0.55_wp,  0.71_wp,       &
    1178                         0.81_wp,  0.90_wp,  0.94_wp,  1.0_wp /)
    1179 
    1180        ENDIF
    1181 
    1182        DO  k = 1, 8
    1183           IF ( collected_r(k) <= mean_rm )  i = k
    1184        ENDDO
    1185 
    1186        DO  k = 1, 18
    1187           IF ( collector_r(k) <= rm )  j = k
    1188        ENDDO
    1189 
    1190        IF ( rm < 10.0_wp )  THEN
    1191           e = 0.0_wp
    1192        ELSEIF ( mean_rm < 2.0_wp )  THEN
    1193           e = 0.001_wp
    1194        ELSEIF ( mean_rm >= 25.0_wp )  THEN
    1195           IF( j <= 2 )  e = 0.0_wp
    1196           IF( j == 3 )  e = 0.47_wp
    1197           IF( j == 4 )  e = 0.8_wp
    1198           IF( j == 5 )  e = 0.9_wp
    1199           IF( j >=6  )  e = 1.0_wp
    1200        ELSEIF ( rm >= 3000.0_wp )  THEN
    1201           IF( i == 1 )  e = 0.02_wp
    1202           IF( i == 2 )  e = 0.16_wp
    1203           IF( i == 3 )  e = 0.33_wp
    1204           IF( i == 4 )  e = 0.55_wp
    1205           IF( i == 5 )  e = 0.71_wp
    1206           IF( i == 6 )  e = 0.81_wp
    1207           IF( i == 7 )  e = 0.90_wp
    1208           IF( i >= 8 )  e = 0.94_wp
    1209        ELSE
    1210           x  = mean_rm - collected_r(i)
    1211           y  = rm - collector_r(j)
    1212           dx = collected_r(i+1) - collected_r(i)
    1213           dy = collector_r(j+1) - collector_r(j)
    1214           aa = x**2 + y**2
    1215           bb = ( dx - x )**2 + y**2
    1216           cc = x**2 + ( dy - y )**2
    1217           dd = ( dx - x )**2 + ( dy - y )**2
    1218           gg = aa + bb + cc + dd
    1219 
    1220           e = ( (gg-aa)*ef(i,j) + (gg-bb)*ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
    1221                 (gg-dd)*ef(i+1,j+1) ) / (3.0_wp * gg)
    1222        ENDIF
    1223 
    1224     END SUBROUTINE collision_efficiency_rogers
    1225 
    1226997 END MODULE lpm_collision_kernels_mod
  • palm/trunk/SOURCE/lpm_data_output_particles.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Tails removed. Unused variables removed.
    2222!
    2323! Former revisions:
     
    5858
    5959    USE control_parameters,                                                    &
    60         ONLY:  prt_time_count, simulated_time
     60        ONLY:  simulated_time
    6161
    6262    USE cpulog,                                                                &
     
    6868    USE kinds
    6969
    70     USE netcdf_interface,                                                      &
    71         ONLY:  nc_stat
    72 
    7370    USE particle_attributes,                                                   &
    74         ONLY:  grid_particles, maximum_number_of_particles,                    &
    75                maximum_number_of_tailpoints, maximum_number_of_tails,          &
    76                number_of_particles, number_of_tails, particles,                &
    77                particle_tail_coordinates, prt_count
     71        ONLY:  grid_particles, number_of_particles,  particles, prt_count
    7872
    7973    IMPLICIT NONE
     
    10397       ENDDO
    10498    ENDDO
    105 !
    106 !-- particle tails currently not available
    107 !    WRITE ( 85 )  maximum_number_of_tailpoints, maximum_number_of_tails, &
    108 !                  number_of_tails
    109 !    IF ( maximum_number_of_tails > 0 )  THEN
    110 !       WRITE ( 85 )  particle_tail_coordinates, prt_time_count
    111 !    ENDIF
    11299
    113100    CALL close_file( 85 )
  • palm/trunk/SOURCE/lpm_droplet_collision.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Integration of a new collision algortithm based on Shima et al. (2009) and
     22! Soelch and Kaercher (2010) called all_or_nothing. The previous implemented
     23! collision algorithm is called average_impact. Moreover, both algorithms are
     24! now positive definit due to their construction, i.e., no negative weighting
     25! factors should occur.
    2226!
    2327! Former revisions:
     
    6367!> Calculates change in droplet radius by collision. Droplet collision is
    6468!> calculated for each grid box seperately. Collision is parameterized by
    65 !> using collision kernels. Three different kernels are available:
    66 !> PALM kernel: Kernel is approximated using a method from Rogers and
    67 !>              Yau (1989, A Short Course in Cloud Physics, Pergamon Press).
    68 !>              All droplets smaller than the treated one are represented by
    69 !>              one droplet with mean features. Collision efficiencies are taken
    70 !>              from the respective table in Rogers and Yau.
     69!> using collision kernels. Two different kernels are available:
    7170!> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which
    7271!>              considers collision due to pure gravitational effects.
     
    8584
    8685    USE arrays_3d,                                                             &
    87         ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
     86        ONLY:  diss, ql_v, ql_vp
    8887
    8988    USE cloud_parameters,                                                      &
    90         ONLY:  effective_coll_efficiency
     89        ONLY:  rho_l
    9190
    9291    USE constants,                                                             &
     
    9493
    9594    USE control_parameters,                                                    &
    96         ONLY:  dt_3d, message_string, u_gtrans, v_gtrans, dz
     95        ONLY:  dt_3d, message_string, dz
    9796
    9897    USE cpulog,                                                                &
     
    10099
    101100    USE grid_variables,                                                        &
    102         ONLY:  ddx, dx, ddy, dy
    103 
    104     USE indices,                                                               &
    105         ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     101        ONLY:  dx, dy
    106102
    107103    USE kinds
    108104
    109105    USE lpm_collision_kernels_mod,                                             &
    110         ONLY:  ckernel, collision_efficiency_rogers, recalculate_kernel
     106        ONLY:  ckernel, recalculate_kernel
    111107
    112108    USE particle_attributes,                                                   &
    113         ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
    114                palm_kernel, particles, particle_type,                          &
    115                prt_count, use_kernel_tables, wang_kernel
     109        ONLY:  all_or_nothing, average_impact, dissipation_classes,            &
     110               hall_kernel, iran_part, number_of_particles, particles,         &
     111               particle_type, prt_count, use_kernel_tables, wang_kernel
     112
     113    USE random_function_mod,                                                   &
     114        ONLY:  random_function
    116115
    117116    USE pegrid
     
    121120    INTEGER(iwp) ::  eclass   !<
    122121    INTEGER(iwp) ::  i        !<
    123     INTEGER(iwp) ::  ii       !<
    124     INTEGER(iwp) ::  inc      !<
    125     INTEGER(iwp) ::  is       !<
    126122    INTEGER(iwp) ::  j        !<
    127     INTEGER(iwp) ::  jj       !<
    128     INTEGER(iwp) ::  js       !<
    129123    INTEGER(iwp) ::  k        !<
    130     INTEGER(iwp) ::  kk       !<
    131124    INTEGER(iwp) ::  n        !<
    132     INTEGER(iwp) ::  pse      !<
    133     INTEGER(iwp) ::  psi      !<
     125    INTEGER(iwp) ::  m        !<
    134126    INTEGER(iwp) ::  rclass_l !<
    135127    INTEGER(iwp) ::  rclass_s !<
    136128
    137     INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !<
    138 
    139     LOGICAL, SAVE ::  first_flag = .TRUE. !<
    140 
    141     TYPE(particle_type) :: tmp_particle   !<
    142 
    143     REAL(wp) ::  aa       !<
    144     REAL(wp) ::  auxn     !< temporary variables
    145     REAL(wp) ::  auxs     !< temporary variables
    146     REAL(wp) ::  bb       !<
    147     REAL(wp) ::  cc       !<
    148     REAL(wp) ::  dd       !<
    149     REAL(wp) ::  ddV      !<
    150     REAL(wp) ::  delta_r  !<
    151     REAL(wp) ::  delta_v  !<
    152     REAL(wp) ::  epsilon  !<
    153     REAL(wp) ::  gg       !<
    154     REAL(wp) ::  mean_r   !<
    155     REAL(wp) ::  ql_int   !<
    156     REAL(wp) ::  ql_int_l !<
    157     REAL(wp) ::  ql_int_u !<
    158     REAL(wp) ::  r3       !<
    159     REAL(wp) ::  sl_r3    !<
    160     REAL(wp) ::  sl_r4    !<
    161     REAL(wp) ::  sum1     !<
    162     REAL(wp) ::  sum2     !<
    163     REAL(wp) ::  sum3     !<
    164     REAL(wp) ::  u_int    !<
    165     REAL(wp) ::  u_int_l  !<
    166     REAL(wp) ::  u_int_u  !<
    167     REAL(wp) ::  v_int    !<
    168     REAL(wp) ::  v_int_l  !<
    169     REAL(wp) ::  v_int_u  !<
    170     REAL(wp) ::  w_int    !<
    171     REAL(wp) ::  w_int_l  !<
    172     REAL(wp) ::  w_int_u  !<
    173     REAL(wp) ::  x        !<
    174     REAL(wp) ::  y        !<
    175 
    176     REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad    !<
    177     REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !<
    178 
    179     REAL, DIMENSION(prt_count(k,j,i))    :: ck
    180     REAL, DIMENSION(prt_count(k,j,i))    :: r3v
    181     REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
    182     REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
     129    REAL(wp) ::  collection_probability  !< probability for collection
     130    REAL(wp) ::  ddV                     !< inverse grid box volume
     131    REAL(wp) ::  epsilon                 !< dissipation rate
     132    REAL(wp) ::  factor_volume_to_mass   !< 4.0 / 3.0 * pi * rho_l
     133    REAL(wp) ::  xm                      !< mean mass of droplet m
     134    REAL(wp) ::  xn                      !< mean mass of droplet n
     135
     136    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight  !< weighting factor
     137    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass    !< total mass of super droplet
    183138
    184139    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
    185140
    186 !
    187 !-- Collision requires at least two particles in the box
    188     IF ( prt_count(k,j,i) > 1 )  THEN
    189 !
    190 !--    First, sort particles within the gridbox by their size,
    191 !--    using Shell's method (see Numerical Recipes)
    192 !--    NOTE: In case of using particle tails, the re-sorting of
    193 !--    ----  tails would have to be included here!
    194        IF ( .NOT. ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    195        use_kernel_tables ) )  THEN
    196           psi = 0
    197           inc = 1
    198           DO WHILE ( inc <= prt_count(k,j,i) )
    199              inc = 3 * inc + 1
    200           ENDDO
    201 
    202           DO WHILE ( inc > 1 )
    203              inc = inc / 3
    204              DO  is = inc+1, prt_count(k,j,i)
    205                 tmp_particle = particles(psi+is)
    206                 js = is
    207                 DO WHILE ( particles(psi+js-inc)%radius >             &
    208                 tmp_particle%radius )
    209                    particles(psi+js) = particles(psi+js-inc)
    210                    js = js - inc
    211                    IF ( js <= inc )  EXIT
    212                 ENDDO
    213                 particles(psi+js) = tmp_particle
    214              ENDDO
    215           ENDDO
    216        ENDIF
    217 
    218        psi = 1
    219        pse = prt_count(k,j,i)
     141    number_of_particles   = prt_count(k,j,i)
     142    factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l
     143    ddV                   = 1 / ( dx * dy * dz )
     144!
     145!-- Collision requires at least one super droplet inside the box
     146    IF ( number_of_particles > 0 )  THEN
    220147
    221148!
    222149!--    Now apply the different kernels
    223        IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    224             use_kernel_tables )  THEN
    225 !
    226 !--       Fast method with pre-calculated efficiencies for
     150       IF ( use_kernel_tables )  THEN
     151!
     152!--       Fast method with pre-calculated collection kernels for
    227153!--       discrete radius- and dissipation-classes.
    228154!--
     
    244170!--       Droplet collision are calculated using collision-coalescence
    245171!--       formulation proposed by Wang (see PALM documentation)
    246 !--       Since new radii after collision are defined by radii of all
    247 !--       droplets before collision, temporary fields for new radii and
    248 !--       weighting factors are needed
    249           ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
    250 
    251           rad    = 0.0_wp
    252           weight = 0.0_wp
    253 
    254           sum1v(1:prt_count(k,j,i)) = 0.0_wp
    255           sum2v(1:prt_count(k,j,i)) = 0.0_wp
    256 
    257           DO  n = 1, prt_count(k,j,i)
    258 
    259              rclass_l = particles(n)%class
    260 !
    261 !--          Mass added due to collisions with smaller droplets
    262              DO  is = n+1, prt_count(k,j,i)
    263                 rclass_s = particles(is)%class
    264                 auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor
    265                 auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor
    266                 IF ( particles(is)%radius <  particles(n)%radius )  THEN
    267                    sum1v(n) =  sum1v(n)  + particles(is)%radius**3 * auxs
    268                    sum2v(is) = sum2v(is) + auxn
    269                 ELSE
    270                    sum2v(n)  = sum2v(n)  + auxs
    271                    sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn
    272                 ENDIF
     172!--       Temporary fields for total mass of super-droplet and weighting factors
     173!--       are allocated.
     174          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
     175
     176          mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
     177                                          particles(1:number_of_particles)%radius**3     * &
     178                                          factor_volume_to_mass
     179          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
     180
     181          IF ( average_impact )  THEN  ! select collision algorithm
     182
     183             DO  n = 1, number_of_particles
     184
     185                rclass_l = particles(n)%class
     186                xn       = mass(n) / weight(n)
     187
     188                DO  m = n, number_of_particles
     189
     190                   rclass_s = particles(m)%class
     191                   xm       = mass(m) / weight(m)
     192
     193                   IF ( xm .LT. xn )  THEN
     194                     
     195!
     196!--                   Particle n collects smaller particle m
     197                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     198                                               weight(n) * ddV * dt_3d
     199
     200                      mass(n)   = mass(n)   + mass(m)   * collection_probability
     201                      weight(m) = weight(m) - weight(m) * collection_probability
     202                      mass(m)   = mass(m)   - mass(m)   * collection_probability
     203                   ELSEIF ( xm .GT. xn )  THEN
     204!
     205!--                   Particle m collects smaller particle n
     206                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     207                                               weight(m) * ddV * dt_3d
     208
     209                      mass(m)   = mass(m)   + mass(n)   * collection_probability
     210                      weight(n) = weight(n) - weight(n) * collection_probability
     211                      mass(n)   = mass(n)   - mass(n)   * collection_probability
     212                   ELSE
     213!
     214!--                   Same-size collections. If n = m, weight is reduced by the
     215!--                   number of possible same-size collections; the total mass
     216!--                   is not changed during same-size collection.
     217!--                   Same-size collections of different
     218!--                   particles ( n /= m ) are treated as same-size collections
     219!--                   of ONE partilce with weight = weight(n) + weight(m) and
     220!--                   mass = mass(n) + mass(m).
     221!--                   Accordingly, each particle loses the same number of
     222!--                   droplets to the other particle, but this has no effect on
     223!--                   total mass mass, since the exchanged droplets have the
     224!--                   same radius.
     225
     226!--                   Note: For m = n this equation is an approximation only
     227!--                   valid for weight >> 1 (which is usually the case). The
     228!--                   approximation is weight(n)-1 = weight(n).
     229                      weight(n) = weight(n) - 0.5_wp * weight(n) *                &
     230                                              ckernel(rclass_l,rclass_s,eclass) * &
     231                                              weight(m) * ddV * dt_3d
     232                      IF ( n .NE. m )  THEN
     233                         weight(m) = weight(m) - 0.5_wp * weight(m) *                &
     234                                                 ckernel(rclass_l,rclass_s,eclass) * &
     235                                                 weight(n) * ddV * dt_3d
     236                      ENDIF
     237                   ENDIF
     238
     239                ENDDO
     240
     241                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     242
    273243             ENDDO
    274           ENDDO
    275           rclass_v = particles(1:prt_count(k,j,i))%class
    276           DO  n = 1, prt_count(k,j,i)
    277             ck(n)       = ckernel(rclass_v(n),rclass_v(n),eclass)
    278           ENDDO
    279           r3v      = particles(1:prt_count(k,j,i))%radius**3
    280           DO  n = 1, prt_count(k,j,i)
    281              sum3 = 0.0_wp
    282              ddV = ddx * ddy / dz
    283 !
    284 !--          Change of the current weighting factor
    285              sum3 = 1 - dt_3d * ddV *                                 &
    286                         ck(n) *                                       &
    287                         ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
    288                     dt_3d * ddV * sum2v(n)
    289              weight(n) = particles(n)%weight_factor * sum3
    290 !
    291 !--          Change of the current droplet radius
    292              rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/&
    293                               sum3 )**0.33333333333333_wp
    294 
    295              ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) &
    296                                          * rad(n)**3
    297 
    298           ENDDO
     244
     245          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
     246
     247             DO  n = 1, number_of_particles
     248
     249                rclass_l = particles(n)%class
     250                xn       = mass(n) / weight(n) ! mean mass of droplet n
     251
     252                DO  m = n, number_of_particles
     253
     254                   rclass_s = particles(m)%class
     255                   xm = mass(m) / weight(m) ! mean mass of droplet m
     256
     257                   IF ( weight(n) .LT. weight(m) )  THEN
     258!
     259!--                   Particle n collects weight(n) droplets of particle m 
     260                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     261                                               weight(m) * ddV * dt_3d
     262
     263                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     264                         mass(n)   = mass(n)   + weight(n) * xm
     265                         weight(m) = weight(m) - weight(n)
     266                         mass(m)   = mass(m)   - weight(n) * xm
     267                      ENDIF
     268
     269                   ELSEIF ( weight(m) .LT. weight(n) )  THEN
     270!
     271!--                   Particle m collects weight(m) droplets of particle n 
     272                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     273                                               weight(n) * ddV * dt_3d
     274
     275                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     276                         mass(m)   = mass(m)   + weight(m) * xn
     277                         weight(n) = weight(n) - weight(m)
     278                         mass(n)   = mass(n)   - weight(m) * xn
     279                      ENDIF
     280                   ELSE
     281!
     282!--                   Collisions of particles of the same weighting factor.
     283!--                   Particle n collects 1/2 weight(n) droplets of particle m,
     284!--                   particle m collects 1/2 weight(m) droplets of particle n.
     285!--                   The total mass mass changes accordingly.
     286!--                   If n = m, the first half of the droplets coalesces with the
     287!--                   second half of the droplets; mass is unchanged because
     288!--                   xm = xn for n = m.
     289
     290!--                   Note: For m = n this equation is an approximation only
     291!--                   valid for weight >> 1 (which is usually the case). The
     292!--                   approximation is weight(n)-1 = weight(n).
     293                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     294                                               weight(n) * ddV * dt_3d
     295
     296                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     297                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
     298                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
     299                         weight(n) = weight(n) - 0.5_wp * weight(m)
     300                         weight(m) = weight(n)
     301                      ENDIF
     302                   ENDIF
     303
     304                ENDDO
     305
     306                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     307
     308             ENDDO
     309
     310          ENDIF
     311
     312
     313
     314
    299315          IF ( ANY(weight < 0.0_wp) )  THEN
    300316                WRITE( message_string, * ) 'negative weighting'
     
    303319          ENDIF
    304320
    305           particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
    306           particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
    307 
    308           DEALLOCATE(rad, weight)
    309 
    310        ELSEIF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    311                 .NOT. use_kernel_tables )  THEN
    312 !
    313 !--       Collision efficiencies are calculated for every new
     321          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
     322                                                      ( weight(1:number_of_particles) &
     323                                                        * factor_volume_to_mass       &
     324                                                      )                               &
     325                                                    )**0.33333333333333_wp
     326
     327          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
     328
     329          DEALLOCATE(weight, mass)
     330
     331       ELSEIF ( .NOT. use_kernel_tables )  THEN
     332!
     333!--       Collection kernels are calculated for every new
    314334!--       grid box. First, allocate memory for kernel table.
    315335!--       Third dimension is 1, because table is re-calculated for
    316336!--       every new dissipation value.
    317           ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) )
    318 !
    319 !--       Now calculate collision efficiencies for this box
     337          ALLOCATE( ckernel(1:number_of_particles,1:number_of_particles,1:1) )
     338!
     339!--       Now calculate collection kernel for this box. Note that
     340!--       the kernel is based on the previous time step
    320341          CALL recalculate_kernel( i, j, k )
    321 
    322342!
    323343!--       Droplet collision are calculated using collision-coalescence
    324344!--       formulation proposed by Wang (see PALM documentation)
    325 !--       Since new radii after collision are defined by radii of all
    326 !--       droplets before collision, temporary fields for new radii and
    327 !--       weighting factors are needed
    328           ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
    329 
    330           rad = 0.0_wp
    331           weight = 0.0_wp
    332 
    333           DO  n = psi, pse, 1
    334 
    335              sum1 = 0.0_wp
    336              sum2 = 0.0_wp
    337              sum3 = 0.0_wp
    338 !
    339 !--          Mass added due to collisions with smaller droplets
    340              DO  is = psi, n-1
    341                 sum1 = sum1 + ( particles(is)%radius**3 *            &
    342                                 ckernel(n,is,1) *  &
    343                                 particles(is)%weight_factor )
     345!--       Temporary fields for total mass of super-droplet and weighting factors
     346!--       are allocated.
     347          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
     348
     349          mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
     350                                        particles(1:number_of_particles)%radius**3     * &
     351                                        factor_volume_to_mass
     352
     353          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
     354
     355          IF ( average_impact )  THEN  ! select collision algorithm
     356
     357             DO  n = 1, number_of_particles
     358
     359                xn = mass(n) / weight(n) ! mean mass of droplet n
     360
     361                DO  m = n, number_of_particles
     362
     363                   xm = mass(m) / weight(m) !mean mass of droplet m
     364
     365                   IF ( xm .LT. xn )  THEN
     366!
     367!--                   Particle n collects smaller particle m
     368                      collection_probability = ckernel(n,m,1) * weight(n) *    &
     369                                               ddV * dt_3d
     370
     371                      mass(n)   = mass(n)   + mass(m)   * collection_probability
     372                      weight(m) = weight(m) - weight(m) * collection_probability
     373                      mass(m)   = mass(m)   - mass(m)   * collection_probability
     374                   ELSEIF ( xm .GT. xn )  THEN
     375!
     376!--                   Particle m collects smaller particle n
     377                      collection_probability = ckernel(n,m,1) * weight(m) *    &
     378                                               ddV * dt_3d
     379
     380                      mass(m)   = mass(m)   + mass(n)   * collection_probability
     381                      weight(n) = weight(n) - weight(n) * collection_probability
     382                      mass(n)   = mass(n)   - mass(n)   * collection_probability
     383                   ELSE
     384!
     385!--                   Same-size collections. If n = m, weight is reduced by the
     386!--                   number of possible same-size collections; the total mass
     387!--                   mass is not changed during same-size collection.
     388!--                   Same-size collections of different
     389!--                   particles ( n /= m ) are treated as same-size collections
     390!--                   of ONE partilce with weight = weight(n) + weight(m) and
     391!--                   mass = mass(n) + mass(m).
     392!--                   Accordingly, each particle loses the same number of
     393!--                   droplets to the other particle, but this has no effect on
     394!--                   total mass mass, since the exchanged droplets have the
     395!--                   same radius.
     396!--
     397!--                   Note: For m = n this equation is an approximation only
     398!--                   valid for weight >> 1 (which is usually the case). The
     399!--                   approximation is weight(n)-1 = weight(n).
     400                      weight(n) = weight(n) - 0.5_wp * weight(n) *             &
     401                                              ckernel(n,m,1) * weight(m) *     &
     402                                              ddV * dt_3d
     403                      IF ( n .NE. m )  THEN
     404                         weight(m) = weight(m) - 0.5_wp * weight(m) *          &
     405                                                 ckernel(n,m,1) * weight(n) *  &
     406                                                 ddV * dt_3d
     407                      ENDIF
     408                   ENDIF
     409
     410
     411                ENDDO
     412
     413                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     414
    344415             ENDDO
    345 !
    346 !--          Rate of collisions with larger droplets
    347              DO  is = n+1, pse
    348                 sum2 = sum2 + ( ckernel(n,is,1) *  &
    349                                 particles(is)%weight_factor )
     416
     417          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
     418
     419             DO  n = 1, number_of_particles
     420
     421                xn = mass(n) / weight(n) ! mean mass of droplet n
     422
     423                DO  m = n, number_of_particles
     424
     425                   xm = mass(m) / weight(m) !mean mass of droplet m
     426
     427                   IF ( weight(n) .LT. weight(m) )  THEN
     428!
     429!--                   Particle n collects smaller particle m
     430                      collection_probability = ckernel(n,m,1) * weight(m) *    &
     431                                               ddV * dt_3d
     432
     433                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     434                         mass(n) = mass(n) + weight(n) * xm 
     435                         weight(m)    = weight(m)    - weight(n)
     436                         mass(m) = mass(m) - weight(n) * xm
     437                      ENDIF
     438
     439                   ELSEIF ( weight(m) .LT. weight(n) )  THEN
     440!
     441!--                   Particle m collects smaller particle n
     442                      collection_probability = ckernel(n,m,1) * weight(n) *    &
     443                                               ddV * dt_3d
     444
     445                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     446                         mass(m) = mass(m) + weight(m) * xn
     447                         weight(n)    = weight(n)    - weight(m)
     448                         mass(n) = mass(n) - weight(m) * xn
     449                      ENDIF
     450                   ELSE
     451!
     452!--                   Collisions of particles of the same weighting factor.
     453!--                   Particle n collects 1/2 weight(n) droplets of particle m,
     454!--                   particle m collects 1/2 weight(m) droplets of particle n.
     455!--                   The total mass mass changes accordingly.
     456!--                   If n = m, the first half of the droplets coalesces with the
     457!--                   second half of the droplets; mass is unchanged because
     458!--                   xm = xn for n = m.
     459!--
     460!--                   Note: For m = n this equation is an approximation only
     461!--                   valid for weight >> 1 (which is usually the case). The
     462!--                   approximation is weight(n)-1 = weight(n).
     463                      collection_probability = ckernel(n,m,1) * weight(n) *    &
     464                                               ddV * dt_3d
     465
     466                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     467                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
     468                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
     469                         weight(n) = weight(n) - 0.5_wp * weight(m)
     470                         weight(m) = weight(n)
     471                      ENDIF
     472                   ENDIF
     473
     474
     475                ENDDO
     476
     477                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     478
    350479             ENDDO
    351480
    352              r3 = particles(n)%radius**3
    353              ddV = ddx * ddy / dz
    354              is = 1
    355 !
    356 !--                   Change of the current weighting factor
    357              sum3 = 1 - dt_3d * ddV *                                 &
    358                         ckernel(n,n,1) *           &
    359                         ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
    360                     dt_3d * ddV * sum2
    361              weight(n-is+1) = particles(n)%weight_factor * sum3
    362 !
    363 !--                   Change of the current droplet radius
    364              rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
    365                               sum3 )**0.33333333333333_wp
    366 
    367              IF ( weight(n-is+1) < 0.0_wp )  THEN
    368                 WRITE( message_string, * ) 'negative weighting',      &
    369                                            'factor: ', weight(n-is+1)
    370                 CALL message( 'lpm_droplet_collision', 'PA0037',      &
     481          ENDIF
     482
     483          IF ( ANY(weight < 0.0_wp) )  THEN
     484                WRITE( message_string, * ) 'negative weighting'
     485                CALL message( 'lpm_droplet_collision', 'PA0028',      &
    371486                               2, 2, -1, 6, 1 )
    372              ENDIF
    373 
    374              ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
    375                                          * rad(n-is+1)**3
    376 
    377           ENDDO
    378 
    379           particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
    380           particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
    381 
    382           DEALLOCATE( rad, weight, ckernel )
    383 
    384        ELSEIF ( palm_kernel )  THEN
    385 !
    386 !--       PALM collision kernel
    387 !
    388 !--       Calculate the mean radius of all those particles which
    389 !--       are of smaller size than the current particle and
    390 !--       use this radius for calculating the collision efficiency
    391           DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
    392 
    393              sl_r3 = 0.0_wp
    394              sl_r4 = 0.0_wp
    395 
    396              DO is = n-1, psi, -1
    397                 IF ( particles(is)%radius < particles(n)%radius )  &
    398                 THEN
    399                    sl_r3 = sl_r3 + particles(is)%weight_factor     &
    400                                    * particles(is)%radius**3
    401                    sl_r4 = sl_r4 + particles(is)%weight_factor     &
    402                                    * particles(is)%radius**4
    403                 ENDIF
    404              ENDDO
    405 
    406              IF ( ( sl_r3 ) > 0.0_wp )  THEN
    407                 mean_r = ( sl_r4 ) / ( sl_r3 )
    408 
    409                 CALL collision_efficiency_rogers( mean_r,             &
    410                                            particles(n)%radius,       &
    411                                            effective_coll_efficiency )
    412 
    413              ELSE
    414                 effective_coll_efficiency = 0.0_wp
    415              ENDIF
    416 
    417              IF ( effective_coll_efficiency > 1.0_wp  .OR.            &
    418                   effective_coll_efficiency < 0.0_wp )                &
    419              THEN
    420                 WRITE( message_string, * )  'collision_efficien' , &
    421                           'cy out of range:' ,effective_coll_efficiency
    422                 CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
    423                               2, -1, 6, 1 )
    424              ENDIF
    425 
    426 !
    427 !--          Interpolation of liquid water content
    428              ii = particles(n)%x * ddx
    429              jj = particles(n)%y * ddy
    430              kk = ( particles(n)%z + 0.5_wp * dz ) / dz
    431 
    432              x  = particles(n)%x - ii * dx
    433              y  = particles(n)%y - jj * dy
    434              aa = x**2          + y**2
    435              bb = ( dx - x )**2 + y**2
    436              cc = x**2          + ( dy - y )**2
    437              dd = ( dx - x )**2 + ( dy - y )**2
    438              gg = aa + bb + cc + dd
    439 
    440              ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
    441                                                   ql(kk,jj,ii+1)   &
    442                         + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
    443                                                   ql(kk,jj+1,ii+1) &
    444                         ) / ( 3.0_wp * gg )
    445 
    446              ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
    447                                                 ql(kk+1,jj,ii+1)   &
    448                         + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
    449                                                 ql(kk+1,jj+1,ii+1) &
    450                         ) / ( 3.0_wp * gg )
    451 
    452              ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
    453                                  ( ql_int_u - ql_int_l )
    454 
    455 !
    456 !--          Interpolate u velocity-component
    457              ii = ( particles(n)%x + 0.5_wp * dx ) * ddx
    458              jj =   particles(n)%y * ddy
    459              kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant
    460 
    461              IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1
    462 
    463              x  = particles(n)%x + ( 0.5_wp - ii ) * dx
    464              y  = particles(n)%y - jj * dy
    465              aa = x**2          + y**2
    466              bb = ( dx - x )**2 + y**2
    467              cc = x**2          + ( dy - y )**2
    468              dd = ( dx - x )**2 + ( dy - y )**2
    469              gg = aa + bb + cc + dd
    470 
    471              u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
    472                                                   u(kk,jj,ii+1)    &
    473                        + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
    474                                                   u(kk,jj+1,ii+1)  &
    475                        ) / ( 3.0_wp * gg ) - u_gtrans
    476              IF ( kk+1 == nzt+1 )  THEN
    477                 u_int = u_int_l
    478              ELSE
    479                 u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
    480                                                  u(kk+1,jj,ii+1)   &
    481                           + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
    482                                                  u(kk+1,jj+1,ii+1) &
    483                           ) / ( 3.0_wp * gg ) - u_gtrans
    484                 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
    485                                   * ( u_int_u - u_int_l )
    486              ENDIF
    487 
    488 !
    489 !--          Same procedure for interpolation of the v velocity-component
    490 !--          (adopt index k from u velocity-component)
    491              ii =   particles(n)%x * ddx
    492              jj = ( particles(n)%y + 0.5_wp * dy ) * ddy
    493 
    494              x  = particles(n)%x - ii * dx
    495              y  = particles(n)%y + ( 0.5_wp - jj ) * dy
    496              aa = x**2          + y**2
    497              bb = ( dx - x )**2 + y**2
    498              cc = x**2          + ( dy - y )**2
    499              dd = ( dx - x )**2 + ( dy - y )**2
    500              gg = aa + bb + cc + dd
    501 
    502              v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
    503                                                    v(kk,jj,ii+1)   &
    504                        + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
    505                                                    v(kk,jj+1,ii+1) &
    506                        ) / ( 3.0_wp * gg ) - v_gtrans
    507              IF ( kk+1 == nzt+1 )  THEN
    508                 v_int = v_int_l
    509              ELSE
    510                 v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
    511                                                    v(kk+1,jj,ii+1) &
    512                           + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
    513                                                  v(kk+1,jj+1,ii+1) &
    514                           ) / ( 3.0_wp * gg ) - v_gtrans
    515                 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
    516                                   * ( v_int_u - v_int_l )
    517              ENDIF
    518 
    519 !
    520 !--          Same procedure for interpolation of the w velocity-component
    521 !--          (adopt index i from v velocity-component)
    522              jj = particles(n)%y * ddy
    523              kk = particles(n)%z / dz
    524 
    525              x  = particles(n)%x - ii * dx
    526              y  = particles(n)%y - jj * dy
    527              aa = x**2          + y**2
    528              bb = ( dx - x )**2 + y**2
    529              cc = x**2          + ( dy - y )**2
    530              dd = ( dx - x )**2 + ( dy - y )**2
    531              gg = aa + bb + cc + dd
    532 
    533              w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
    534                                                      w(kk,jj,ii+1) &
    535                        + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
    536                                                    w(kk,jj+1,ii+1) &
    537                        ) / ( 3.0_wp * gg )
    538              IF ( kk+1 == nzt+1 )  THEN
    539                 w_int = w_int_l
    540              ELSE
    541                 w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
    542                                                    w(kk+1,jj,ii+1) &
    543                           + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
    544                                                  w(kk+1,jj+1,ii+1) &
    545                           ) / ( 3.0_wp * gg )
    546                 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
    547                                   * ( w_int_u - w_int_l )
    548              ENDIF
    549 
    550 !
    551 !--          Change in radius due to collision
    552              delta_r = effective_coll_efficiency / 3.0_wp          &
    553                        * pi * sl_r3 * ddx * ddy / dz               &
    554                        * SQRT( ( u_int - particles(n)%speed_x )**2 &
    555                              + ( v_int - particles(n)%speed_y )**2 &
    556                              + ( w_int - particles(n)%speed_z )**2 &
    557                              ) * dt_3d
    558 !
    559 !--          Change in volume due to collision
    560              delta_v = particles(n)%weight_factor                  &
    561                        * ( ( particles(n)%radius + delta_r )**3    &
    562                            - particles(n)%radius**3 )
    563 
    564 !
    565 !--          Check if collected particles provide enough LWC for
    566 !--          volume change of collector particle
    567              IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
    568 
    569                 delta_r = ( ( sl_r3/particles(n)%weight_factor )               &
    570                             + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp )    &
    571                             - particles(n)%radius
    572 
    573                 DO  is = n-1, psi, -1
    574                    IF ( particles(is)%radius < particles(n)%radius )  THEN
    575                       particles(is)%weight_factor = 0.0_wp
    576                       particles(is)%particle_mask = .FALSE.
    577                       deleted_particles = deleted_particles + 1
    578                    ENDIF
    579                 ENDDO
    580 
    581              ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
    582 
    583                 DO  is = n-1, psi, -1
    584                    IF ( particles(is)%radius < particles(n)%radius &
    585                         .AND.  sl_r3 > 0.0_wp )  THEN
    586                       particles(is)%weight_factor =                &
    587                                    ( ( particles(is)%weight_factor &
    588                                    * ( particles(is)%radius**3 ) ) &
    589                                    - ( delta_v                     &
    590                                    * particles(is)%weight_factor   &
    591                                    * ( particles(is)%radius**3 )   &
    592                                    / sl_r3 ) )                     &
    593                                    / ( particles(is)%radius**3 )
    594 
    595                       IF ( particles(is)%weight_factor < 0.0_wp )  THEN
    596                          WRITE( message_string, * ) 'negative ',   &
    597                                          'weighting factor: ',     &
    598                                          particles(is)%weight_factor
    599                          CALL message( 'lpm_droplet_collision', &
    600                                        'PA0039',                &
    601                                        2, 2, -1, 6, 1 )
    602                       ENDIF
    603                    ENDIF
    604                 ENDDO
    605 
    606              ENDIF
    607 
    608              particles(n)%radius = particles(n)%radius + delta_r
    609              ql_vp(k,j,i) = ql_vp(k,j,i) +               &
    610                             particles(n)%weight_factor * &
    611                             ( particles(n)%radius**3 )
    612           ENDDO
    613 
    614           ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    615                                         * particles(psi)%radius**3
    616 
    617        ENDIF  ! collision kernel
    618 
    619     ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    620 
    621        psi = 1
    622 
    623 !
    624 !--    Calculate change of weighting factor due to self collision
    625        IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    626             use_kernel_tables )  THEN
    627 
    628           IF ( wang_kernel )  THEN
    629              eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
    630                            dissipation_classes ) + 1
    631              epsilon = diss(k,j,i)
    632           ELSE
    633              epsilon = 0.0_wp
    634487          ENDIF
    635           IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
    636              eclass = 0   ! Hall kernel is used
    637           ELSE
    638              eclass = MIN( dissipation_classes, eclass )
    639           ENDIF
    640 
    641           ddV = ddx * ddy / dz
    642           rclass_l = particles(psi)%class
    643           sum3 = 1 - dt_3d * ddV *                                 &
    644                      ( ckernel(rclass_l,rclass_l,eclass) *         &
    645                        ( particles(psi)%weight_factor-1 ) * 0.5_wp )
    646 
    647           particles(psi)%radius = ( particles(psi)%radius**3 / &
    648                                     sum3 )**0.33333333333333_wp
    649           particles(psi)%weight_factor = particles(psi)%weight_factor &
    650                                          * sum3
    651 
    652        ELSE IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    653                   .NOT. use_kernel_tables )  THEN
    654 !
    655 !--       Collision efficiencies are calculated for every new
    656 !--       grid box. First, allocate memory for kernel table.
    657 !--       Third dimension is 1, because table is re-calculated for
    658 !--       every new dissipation value.
    659           ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
    660 !
    661 !--       Now calculate collision efficiencies for this box
    662           CALL recalculate_kernel( i, j, k )
    663 
    664           ddV = ddx * ddy / dz
    665           sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
    666                        ( particles(psi)%weight_factor - 1 ) * 0.5_wp )
    667 
    668           particles(psi)%radius = ( particles(psi)%radius**3 / &
    669                                     sum3 )**0.33333333333333_wp
    670           particles(psi)%weight_factor = particles(psi)%weight_factor &
    671                                          * sum3
    672 
    673           DEALLOCATE( ckernel )
    674        ENDIF
    675 
    676       ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
    677                        particles(psi)%radius**3
     488
     489          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
     490                                                      ( weight(1:number_of_particles) &
     491                                                        * factor_volume_to_mass       &
     492                                                      )                               &
     493                                                    )**0.33333333333333_wp
     494
     495          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
     496
     497          DEALLOCATE( weight, mass, ckernel )
     498
     499       ENDIF
     500 
    678501    ENDIF
    679 
    680 !
    681 !-- Check if condensation of LWC was conserved during collision process
     502   
     503
     504!
     505!-- Check if LWC is conserved during collision process
    682506    IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
    683        IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.             &
     507       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.                      &
    684508            ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp )  THEN
    685           WRITE( message_string, * ) 'LWC is not conserved during',&
    686                                      ' collision! ',               &
    687                                      'LWC after condensation: ',   &
    688                                      ql_v(k,j,i),                  &
    689                                      ' LWC after collision: ',     &
    690                                      ql_vp(k,j,i)
    691           CALL message( 'lpm_droplet_collision', 'PA0040',         &
    692                         2, 2, -1, 6, 1 )
     509          WRITE( message_string, * ) ' LWC is not conserved during',           &
     510                                     ' collision! ',                           &
     511                                     ' LWC after condensation: ', ql_v(k,j,i), &
     512                                     ' LWC after collision: ', ql_vp(k,j,i)
     513          CALL message( 'lpm_droplet_collision', 'PA0040', 2, 2, -1, 6, 1 )
    693514       ENDIF
    694515    ENDIF
  • palm/trunk/SOURCE/lpm_droplet_condensation.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Unused variables removed.
    2222!
    2323! Former revisions:
     
    8787
    8888    USE cloud_parameters,                                                      &
    89         ONLY:  bfactor, curvature_solution_effects, diff_coeff_l,              &
    90                eps_ros, l_d_rv, l_v, rho_l,  r_v, thermal_conductivity_l
     89        ONLY:  bfactor, curvature_solution_effects, eps_ros, l_d_rv, l_v,      &
     90               rho_l,  r_v
    9191
    9292    USE constants,                                                             &
     
    9494
    9595    USE control_parameters,                                                    &
    96         ONLY:  atmos_ocean_sign, dt_3d, dz, message_string,                    &
    97                molecular_viscosity, rho_surface
     96        ONLY:  dt_3d, dz, message_string, molecular_viscosity, rho_surface
     97
    9898    USE cpulog,                                                                &
    9999        ONLY:  cpu_log, log_point_s
    100100
    101101    USE grid_variables,                                                        &
    102         ONLY:  dx, ddx, dy, ddy
     102        ONLY:  dx, dy
    103103
    104104    USE lpm_collision_kernels_mod,                                             &
     
    109109    USE particle_attributes,                                                   &
    110110        ONLY:  block_offset, grid_particles, hall_kernel, number_of_particles, &
    111                offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
    112                radius_classes, use_kernel_tables, wang_kernel
     111               particles, radius_classes, use_kernel_tables, wang_kernel
    113112
    114113
  • palm/trunk/SOURCE/lpm_exchange_horiz.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Tails removed. Unused variables removed.
    2222!
    2323! Former revisions:
     
    6767! Description:
    6868! ------------
    69 !> Exchange of particles (and tails) between the subdomains.
     69! Exchange of particles between the subdomains.
    7070!------------------------------------------------------------------------------!
    7171 MODULE lpm_exchange_horiz_mod
     
    9393
    9494    USE particle_attributes,                                                   &
    95         ONLY:  alloc_factor, deleted_particles, deleted_tails, grid_particles, &
    96                ibc_par_lr, ibc_par_ns, maximum_number_of_tails,                &
    97                maximum_number_of_tailpoints, min_nr_particle,                  &
    98                mpi_particle_type, number_of_tails, number_of_particles,        &
    99                offset_ocean_nzt, particles,                                    &
    100                particle_tail_coordinates, particle_type, prt_count,            &
    101                tail_mask, trlp_count_sum,                                      &
     95        ONLY:  alloc_factor, deleted_particles, grid_particles,                &
     96               ibc_par_lr, ibc_par_ns, min_nr_particle,                        &
     97               mpi_particle_type, number_of_particles,                         &
     98               offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
     99               particle_type, prt_count, trlp_count_sum,                       &
    102100               trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
    103101               trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
    104                trsp_count_recv_sum, use_particle_tails, zero_particle
     102               trsp_count_recv_sum, zero_particle
    105103
    106104    USE pegrid
     
    156154    INTEGER(iwp) ::  j                                       !<
    157155    INTEGER(iwp) ::  jp                                      !<
    158     INTEGER(iwp) ::  k                                       !<
    159156    INTEGER(iwp) ::  kp                                      !<
    160157    INTEGER(iwp) ::  n                                       !<
    161     INTEGER(iwp) ::  nn                                      !<
    162     INTEGER(iwp) ::  tlength                                 !<
    163158    INTEGER(iwp) ::  trlp_count                              !<
    164159    INTEGER(iwp) ::  trlp_count_recv                         !<
     
    178173    INTEGER(iwp) ::  trspt_count_recv                        !<
    179174
    180     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trlpt        !<
    181     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trnpt        !<
    182     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trrpt        !<
    183     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trspt        !<
    184 
    185175    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !<
    186176    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !<
     
    196186#if defined( __parallel )
    197187
     188!
     189!-- Exchange between subdomains.
     190!-- As soon as one particle has moved beyond the boundary of the domain, it
     191!-- is included in the relevant transfer arrays and marked for subsequent
     192!-- deletion on this PE.
     193!-- First sweep for crossings in x direction. Find out first the number of
     194!-- particles to be transferred and allocate temporary arrays needed to store
     195!-- them.
     196!-- For a one-dimensional decomposition along y, no transfer is necessary,
     197!-- because the particle remains on the PE, but the particle coordinate has to
     198!-- be adjusted.
    198199    trlp_count  = 0
    199200    trlpt_count = 0
     
    226227                      IF ( i < nxl )  THEN
    227228                         trlp_count = trlp_count + 1
    228                          IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
    229229                      ELSEIF ( i > nxr )  THEN
    230230                         trrp_count = trrp_count + 1
    231                          IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
    232231                      ENDIF
    233232                   ENDIF
     
    247246       trlp = zero_particle
    248247       trrp = zero_particle
    249 
    250        IF ( use_particle_tails )  THEN
    251           ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
    252                     trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
    253           tlength = maximum_number_of_tailpoints * 5
    254        ENDIF
    255248
    256249       trlp_count  = 0
     
    269262             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    270263             DO  n = 1, number_of_particles
    271 
    272                 nn = particles(n)%tail_id
    273264!
    274265!--             Only those particles that have not been marked as 'deleted' may
     
    292283                               particles(n)%origin_x = ( nx + 1 ) * dx + &
    293284                               particles(n)%origin_x
    294                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    295                                   i  = particles(n)%tailpoints
    296                                   particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
    297                                   + particle_tail_coordinates(1:i,1,nn)
    298                                ENDIF
    299285                            ELSE
    300286                               trlp_count = trlp_count + 1
     
    312298                               ENDIF
    313299
    314                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    315                                   trlpt_count = trlpt_count + 1
    316                                   trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
    317                                   trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
    318                                   trlpt(:,1,trlpt_count)
    319                                   tail_mask(nn) = .FALSE.
    320                                   deleted_tails = deleted_tails + 1
    321                                ENDIF
    322300                            ENDIF
    323301
     
    327305                            particles(n)%particle_mask = .FALSE.
    328306                            deleted_particles = deleted_particles + 1
    329                             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    330                                tail_mask(nn) = .FALSE.
    331                                deleted_tails = deleted_tails + 1
    332                             ENDIF
    333307
    334308                         ELSEIF ( ibc_par_lr == 2 )  THEN
     
    348322                         deleted_particles = deleted_particles + 1
    349323
    350                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    351                             trlpt_count = trlpt_count + 1
    352                             trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
    353                             tail_mask(nn) = .FALSE.
    354                             deleted_tails = deleted_tails + 1
    355                          ENDIF
    356324                      ENDIF
    357325
     
    367335                               particles(n)%origin_x = particles(n)%origin_x - &
    368336                               ( nx + 1 ) * dx
    369                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    370                                   i = particles(n)%tailpoints
    371                                   particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
    372                                   + particle_tail_coordinates(1:i,1,nn)
    373                                ENDIF
    374337                            ELSE
    375338                               trrp_count = trrp_count + 1
     
    381344                               deleted_particles = deleted_particles + 1
    382345
    383                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    384                                   trrpt_count = trrpt_count + 1
    385                                   trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
    386                                   trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
    387                                   ( nx + 1 ) * dx
    388                                   tail_mask(nn) = .FALSE.
    389                                   deleted_tails = deleted_tails + 1
    390                                ENDIF
    391346                            ENDIF
    392347
     
    396351                            particles(n)%particle_mask = .FALSE.
    397352                            deleted_particles = deleted_particles + 1
    398                             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    399                                tail_mask(nn) = .FALSE.
    400                                deleted_tails = deleted_tails + 1
    401                             ENDIF
    402353
    403354                         ELSEIF ( ibc_par_lr == 2 )  THEN
     
    417368                         deleted_particles = deleted_particles + 1
    418369
    419                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    420                             trrpt_count = trrpt_count + 1
    421                             trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
    422                             tail_mask(nn) = .FALSE.
    423                             deleted_tails = deleted_tails + 1
    424                          ENDIF
    425370                      ENDIF
    426371
     
    452397       DEALLOCATE(rvrp)
    453398
    454        IF ( use_particle_tails )  THEN
    455 
    456           CALL MPI_SENDRECV( trlpt_count,      1, MPI_INTEGER, pleft,  0, &
    457                              trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
    458                              comm2d, status, ierr )
    459 
    460           IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
    461           THEN
    462              IF ( netcdf_data_format < 3 )  THEN
    463                 message_string = 'maximum_number_of_tails ' //   &
    464                                  'needs to be increased ' //     &
    465                                  '&but this is not allowed wi'// &
    466                                  'th netcdf_data_format < 3'
    467                 CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
    468              ELSE
    469                 CALL lpm_extend_tail_array( trrpt_count_recv )
    470              ENDIF
    471           ENDIF
    472 
    473           CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,      &
    474                              pleft, 1,                                         &
    475                              particle_tail_coordinates(1,1,number_of_tails+1), &
    476                              trrpt_count_recv*tlength, MPI_REAL, pright, 1,    &
    477                              comm2d, status, ierr )
    478 !
    479 !--       Update the tail ids for the transferred particles
    480           nn = number_of_tails
    481           DO  n = number_of_particles+1, number_of_particles+trrp_count_recv
    482              IF ( particles(n)%tail_id /= 0 )  THEN
    483                 nn = nn + 1
    484                 particles(n)%tail_id = nn
    485              ENDIF
    486           ENDDO
    487 
    488        ENDIF
    489 
    490399!
    491400!--    Send right boundary, receive left boundary
     
    503412       IF ( trlp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
    504413
    505        DEALLOCATE(rvlp)
    506 
    507        IF ( use_particle_tails )  THEN
    508 
    509           CALL MPI_SENDRECV( trrpt_count,      1, MPI_INTEGER, pright, 0, &
    510                              trlpt_count_recv, 1, MPI_INTEGER, pleft,  0, &
    511                              comm2d, status, ierr )
    512 
    513           IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
    514           THEN
    515              IF ( netcdf_data_format < 3 )  THEN
    516                 message_string = 'maximum_number_of_tails ' //   &
    517                                  'needs to be increased ' //     &
    518                                  '&but this is not allowed wi'// &
    519                                  'th netcdf_data_format < 3'
    520                 CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
    521              ELSE
    522                 CALL lpm_extend_tail_array( trlpt_count_recv )
    523              ENDIF
    524           ENDIF
    525 
    526           CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,      &
    527                              pright, 1,                                        &
    528                              particle_tail_coordinates(1,1,number_of_tails+1), &
    529                              trlpt_count_recv*tlength, MPI_REAL, pleft, 1,     &
    530                              comm2d, status, ierr )
    531 !
    532 !--       Update the tail ids for the transferred particles
    533           nn = number_of_tails
    534           DO  n = number_of_particles+1, number_of_particles+trlp_count_recv
    535              IF ( particles(n)%tail_id /= 0 )  THEN
    536                 nn = nn + 1
    537                 particles(n)%tail_id = nn
    538              ENDIF
    539           ENDDO
    540 
    541        ENDIF
    542 
    543 !       number_of_particles = number_of_particles + trlp_count_recv
    544 !       number_of_tails     = number_of_tails     + trlpt_count_recv
    545 
    546        IF ( use_particle_tails )  THEN
    547           DEALLOCATE( trlpt, trrpt )
    548        ENDIF
     414       DEALLOCATE( rvlp )
    549415       DEALLOCATE( trlp, trrp )
    550416
     
    589455                      IF ( j < nys )  THEN
    590456                         trsp_count = trsp_count + 1
    591                          IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count + 1
    592457                      ELSEIF ( j > nyn )  THEN
    593458                         trnp_count = trnp_count + 1
    594                          IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count + 1
    595459                      ENDIF
    596460                   ENDIF
     
    609473       trsp = zero_particle
    610474       trnp = zero_particle
    611 
    612        IF ( use_particle_tails )  THEN
    613           ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
    614                     trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
    615           tlength = maximum_number_of_tailpoints * 5
    616        ENDIF
    617475
    618476       trsp_count  = nr_move_south
     
    633491             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    634492             DO  n = 1, number_of_particles
    635 
    636                 nn = particles(n)%tail_id
    637493!
    638494!--             Only those particles that have not been marked as 'deleted' may
     
    655511                               particles(n)%origin_y = ( ny + 1 ) * dy + &
    656512                               particles(n)%origin_y
    657                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    658                                   i = particles(n)%tailpoints
    659                                   particle_tail_coordinates(1:i,2,nn) =        &
    660                                      ( ny+1 ) * dy + particle_tail_coordinates(1:i,2,nn)
    661                                ENDIF
    662513                            ELSE
    663514                               trsp_count = trsp_count + 1
     
    677528                               ENDIF
    678529
    679                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    680                                   trspt_count = trspt_count + 1
    681                                   trspt(:,:,trspt_count) = &
    682                                   particle_tail_coordinates(:,:,nn)
    683                                   trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
    684                                   trspt(:,2,trspt_count)
    685                                   tail_mask(nn) = .FALSE.
    686                                   deleted_tails = deleted_tails + 1
    687                                ENDIF
    688530                            ENDIF
    689531
     
    693535                            particles(n)%particle_mask = .FALSE.
    694536                            deleted_particles = deleted_particles + 1
    695                             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    696                                tail_mask(nn) = .FALSE.
    697                                deleted_tails = deleted_tails + 1
    698                             ENDIF
    699537
    700538                         ELSEIF ( ibc_par_ns == 2 )  THEN
     
    714552                         deleted_particles = deleted_particles + 1
    715553
    716                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    717                             trspt_count = trspt_count + 1
    718                             trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
    719                             tail_mask(nn) = .FALSE.
    720                             deleted_tails = deleted_tails + 1
    721                          ENDIF
    722554                      ENDIF
    723555
     
    733565                               particles(n)%origin_y =                         &
    734566                                  particles(n)%origin_y - ( ny + 1 ) * dy
    735                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    736                                   i = particles(n)%tailpoints
    737                                   particle_tail_coordinates(1:i,2,nn) =        &
    738                                      - (ny+1) * dy + particle_tail_coordinates(1:i,2,nn)
    739                                ENDIF
    740567                            ELSE
    741568                               trnp_count = trnp_count + 1
     
    747574                               particles(n)%particle_mask = .FALSE.
    748575                               deleted_particles = deleted_particles + 1
    749                                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    750                                   trnpt_count = trnpt_count + 1
    751                                   trnpt(:,:,trnpt_count) =                     &
    752                                      particle_tail_coordinates(:,:,nn)
    753                                   trnpt(:,2,trnpt_count) =                     &
    754                                      trnpt(:,2,trnpt_count) - ( ny + 1 ) * dy
    755                                   tail_mask(nn) = .FALSE.
    756                                   deleted_tails = deleted_tails + 1
    757                                ENDIF
    758576                            ENDIF
    759577
     
    763581                            particles(n)%particle_mask = .FALSE.
    764582                            deleted_particles = deleted_particles + 1
    765                             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    766                                tail_mask(nn) = .FALSE.
    767                                deleted_tails = deleted_tails + 1
    768                             ENDIF
    769583
    770584                         ELSEIF ( ibc_par_ns == 2 )  THEN
     
    784598                         deleted_particles = deleted_particles + 1
    785599
    786                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    787                             trnpt_count = trnpt_count + 1
    788                             trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
    789                             tail_mask(nn) = .FALSE.
    790                             deleted_tails = deleted_tails + 1
    791                          ENDIF
    792600                      ENDIF
    793601
     
    819627       DEALLOCATE(rvnp)
    820628
    821        IF ( use_particle_tails )  THEN
    822 
    823           CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
    824                              trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
    825                              comm2d, status, ierr )
    826 
    827           IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
    828           THEN
    829              IF ( netcdf_data_format < 3 )  THEN
    830                 message_string = 'maximum_number_of_tails ' //    &
    831                                  'needs to be increased ' //      &
    832                                  '&but this is not allowed wi' // &
    833                                  'th netcdf_data_format < 3'
    834                 CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
    835              ELSE
    836                 CALL lpm_extend_tail_array( trnpt_count_recv )
    837              ENDIF
    838           ENDIF
    839 
    840           CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,      &
    841                              psouth, 1,                                        &
    842                              particle_tail_coordinates(1,1,number_of_tails+1), &
    843                              trnpt_count_recv*tlength, MPI_REAL, pnorth, 1,    &
    844                              comm2d, status, ierr )
    845 
    846 !
    847 !--       Update the tail ids for the transferred particles
    848           nn = number_of_tails
    849           DO  n = number_of_particles+1, number_of_particles+trnp_count_recv
    850              IF ( particles(n)%tail_id /= 0 )  THEN
    851                 nn = nn + 1
    852                 particles(n)%tail_id = nn
    853              ENDIF
    854           ENDDO
    855 
    856        ENDIF
    857 
    858 !       number_of_particles = number_of_particles + trnp_count_recv
    859 !       number_of_tails     = number_of_tails     + trnpt_count_recv
    860 
    861629!
    862630!--    Send back boundary, receive front boundary
     
    876644       DEALLOCATE(rvsp)
    877645
    878        IF ( use_particle_tails )  THEN
    879 
    880           CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
    881                              trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
    882                              comm2d, status, ierr )
    883 
    884           IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
    885           THEN
    886              IF ( netcdf_data_format < 3 )  THEN
    887                 message_string = 'maximum_number_of_tails ' //   &
    888                                  'needs to be increased ' //     &
    889                                  '&but this is not allowed wi'// &
    890                                  'th NetCDF output switched on'
    891                 CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
    892              ELSE
    893                 CALL lpm_extend_tail_array( trspt_count_recv )
    894              ENDIF
    895           ENDIF
    896 
    897           CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,      &
    898                              pnorth, 1,                                        &
    899                              particle_tail_coordinates(1,1,number_of_tails+1), &
    900                              trspt_count_recv*tlength, MPI_REAL, psouth, 1,    &
    901                              comm2d, status, ierr )
    902 !
    903 !--       Update the tail ids for the transferred particles
    904           nn = number_of_tails
    905           DO  n = number_of_particles+1, number_of_particles+trsp_count_recv
    906              IF ( particles(n)%tail_id /= 0 )  THEN
    907                 nn = nn + 1
    908                 particles(n)%tail_id = nn
    909              ENDIF
    910           ENDDO
    911 
    912        ENDIF
    913 
    914646       number_of_particles = number_of_particles + trsp_count_recv
    915        number_of_tails     = number_of_tails     + trspt_count_recv
    916 
    917        IF ( use_particle_tails )  THEN
    918           DEALLOCATE( trspt, trnpt )
    919        ENDIF
     647
    920648       DEALLOCATE( trsp, trnp )
    921649
     
    928656    DO  n = 1, number_of_particles
    929657
    930        nn = particles(n)%tail_id
    931 
    932658       IF ( particles(n)%x < -0.5_wp * dx )  THEN
    933659
     
    936662!--          Cyclic boundary. Relevant coordinate has to be changed.
    937663             particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
    938              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    939                 i = particles(n)%tailpoints
    940                 particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
    941                                              particle_tail_coordinates(1:i,1,nn)
    942              ENDIF
     664
    943665          ELSEIF ( ibc_par_lr == 1 )  THEN
    944666!
     
    946668             particles(n)%particle_mask = .FALSE.
    947669             deleted_particles = deleted_particles + 1
    948              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    949                 tail_mask(nn) = .FALSE.
    950                 deleted_tails = deleted_tails + 1
    951              ENDIF
     670
    952671          ELSEIF ( ibc_par_lr == 2 )  THEN
    953672!
     
    963682!--          Cyclic boundary. Relevant coordinate has to be changed.
    964683             particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
    965              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    966                 i = particles(n)%tailpoints
    967                 particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
    968                                              particle_tail_coordinates(1:i,1,nn)
    969              ENDIF
     684
    970685          ELSEIF ( ibc_par_lr == 1 )  THEN
    971686!
     
    973688             particles(n)%particle_mask = .FALSE.
    974689             deleted_particles = deleted_particles + 1
    975              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    976                 tail_mask(nn) = .FALSE.
    977                 deleted_tails = deleted_tails + 1
    978              ENDIF
     690
    979691          ELSEIF ( ibc_par_lr == 2 )  THEN
    980692!
     
    992704!--          Cyclic boundary. Relevant coordinate has to be changed.
    993705             particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
    994              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    995                 i = particles(n)%tailpoints
    996                 particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
    997                                              particle_tail_coordinates(1:i,2,nn)
    998              ENDIF
     706
    999707          ELSEIF ( ibc_par_ns == 1 )  THEN
    1000708!
     
    1002710             particles(n)%particle_mask = .FALSE.
    1003711             deleted_particles = deleted_particles + 1
    1004              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    1005                 tail_mask(nn) = .FALSE.
    1006                 deleted_tails = deleted_tails + 1
    1007              ENDIF
     712
    1008713          ELSEIF ( ibc_par_ns == 2 )  THEN
    1009714!
     
    1019724!--          Cyclic boundary. Relevant coordinate has to be changed.
    1020725             particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
    1021              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    1022                 i = particles(n)%tailpoints
    1023                 particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
    1024                                              particle_tail_coordinates(1:i,2,nn)
    1025              ENDIF
     726
    1026727          ELSEIF ( ibc_par_ns == 1 )  THEN
    1027728!
     
    1029730             particles(n)%particle_mask = .FALSE.
    1030731             deleted_particles = deleted_particles + 1
    1031              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    1032                 tail_mask(nn) = .FALSE.
    1033                 deleted_tails = deleted_tails + 1
    1034              ENDIF
     732
    1035733          ELSEIF ( ibc_par_ns == 2 )  THEN
    1036734!
     
    1287985    INTEGER(iwp), INTENT(in)                       ::  j              !<
    1288986    INTEGER(iwp), INTENT(in)                       ::  k              !<
    1289     INTEGER(iwp), INTENT(in), optional             ::  size_in        !<
     987    INTEGER(iwp), INTENT(in), OPTIONAL             ::  size_in        !<
    1290988
    1291989    INTEGER(iwp)                                   :: old_size        !<
     
    1300998       new_size = size_in
    1301999    ELSE
    1302        new_size = old_size * ( 1.0 + alloc_factor / 100.0 )
     1000       new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
    13031001    ENDIF
    13041002
  • palm/trunk/SOURCE/lpm_init.f90

    r1818 r1822  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Unused variables removed.
    2222!
    2323! Former revisions:
     
    113113               dz, initializing_actions, message_string, ocean, simulated_time
    114114
    115     USE dvrp_variables,                                                        &
    116         ONLY:  particle_color, particle_dvrpsize
    117 
    118115    USE grid_variables,                                                        &
    119116        ONLY:  ddx, dx, ddy, dy
     
    134131        ONLY:   alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,        &
    135132                block_offset, block_offset_def, collision_kernel,              &
    136                 density_ratio, dvrp_psize, grid_particles,                     &
     133                density_ratio, grid_particles,                                 &
    137134                initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns,   &
    138135                ibc_par_t, iran_part, log_z_z0,                                &
    139136                max_number_of_particle_groups, maximum_number_of_particles,    &
    140                 maximum_number_of_tailpoints, maximum_number_of_tails,         &
    141                 minimum_tailpoint_distance, min_nr_particle,                   &
    142                 mpi_particle_type, new_tail_id,                                &
    143                 number_of_initial_tails, number_of_particles,                  &
     137                min_nr_particle, mpi_particle_type,                            &
     138                number_of_particles,                                           &
    144139                number_of_particle_groups, number_of_sublayers,                &
    145                 number_of_tails, offset_ocean_nzt, offset_ocean_nzt_m1,        &
     140                offset_ocean_nzt, offset_ocean_nzt_m1,                         &
    146141                particles, particle_advection_start, particle_groups,          &
    147142                particle_groups_type, particles_per_point,                     &
    148                 particle_tail_coordinates,  particle_type, pdx, pdy, pdz,      &
     143                particle_type, pdx, pdy, pdz,                                  &
    149144                prt_count, psb, psl, psn, psr, pss, pst,                       &
    150145                radius, random_start_position, read_particles_from_restartfile,&
    151                 seed_follows_topography, skip_particles_for_tail, sort_count,  &
    152                 tail_mask, total_number_of_particles, total_number_of_tails,   &
    153                 use_particle_tails, use_sgs_for_particles,                     &
     146                seed_follows_topography, sort_count,                           &
     147                total_number_of_particles,                                     &
     148                use_sgs_for_particles,                                         &
    154149                write_particle_statistics, uniform_particles, zero_particle,   &
    155150                z0_av_global
     
    192187
    193188    INTEGER(iwp) ::  i                           !<
    194     INTEGER(iwp) ::  ip                          !<
    195189    INTEGER(iwp) ::  j                           !<
    196     INTEGER(iwp) ::  jp                          !<
    197190    INTEGER(iwp) ::  k                           !<
    198     INTEGER(iwp) ::  kp                          !<
    199     INTEGER(iwp) ::  n                           !<
    200     INTEGER(iwp) ::  nn                          !<
    201191
    202192#if defined( __parallel )
     
    205195    INTEGER(iwp), DIMENSION(3) ::  types         !<
    206196#endif
    207     LOGICAL ::  uniform_particles_l              !<
    208197
    209198    REAL(wp) ::  height_int                      !<
     
    216205!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
    217206!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
    218 #if defined( __twocachelines )
    219     blocklengths(1)  =  7;  blocklengths(2)  =  18;  blocklengths(3)  =   1
    220     displacements(1) =  0;  displacements(2) =  64;  displacements(3) = 128
    221 
    222     types(1) = MPI_REAL                               ! 64 bit words
    223     types(2) = MPI_INTEGER                            ! 32 Bit words
    224     types(3) = MPI_UB
    225 #else
    226207    blocklengths(1)  = 19;  blocklengths(2)  =   6;  blocklengths(3)  =   1
    227208    displacements(1) =  0;  displacements(2) = 152;  displacements(3) = 176
     
    230211    types(2) = MPI_INTEGER
    231212    types(3) = MPI_UB
    232 #endif
    233213    CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
    234214                          mpi_particle_type, ierr )
     
    293273!
    294274!-- Allocate arrays required for calculating particle SGS velocities
    295     IF ( use_sgs_for_particles )  THEN
     275    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    296276       ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    297277                 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     
    429409!--    Allocate particle arrays and set attributes of the initial set of
    430410!--    particles, which can be also periodically released at later times.
    431 !--    Also allocate array for particle tail coordinates, if needed.
    432411       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    433412                 grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
     
    443422!--    occur within restart runs). The reason for this is still not clear
    444423!--    and may be presumably caused by errors in the respective user-interface.
    445 #if defined( __twocachelines )
    446        zero_particle = particle_type( 0.0_wp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp,  &
    447                                       0.0_sp, 0.0_sp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    448                                       0, .FALSE., 0.0_wp, 0.0_wp, 0.0_wp,      &
    449                                       0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp,  &
    450                                       0.0_sp, 0, 0, 0, -1)
    451 #else
    452424       zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    453425                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
     
    455427                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
    456428                                      0, .FALSE., -1)
    457 #endif
     429
    458430       particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
    459 
    460 !
    461 !--    Set the default particle size used for dvrp plots
    462        IF ( dvrp_psize == 9999999.9_wp )  dvrp_psize = 0.2_wp * dx
    463431
    464432!
     
    506474       ENDIF
    507475
    508 !
    509 !--    Check if particles are really uniform in color and radius (dvrp_size)
    510 !--    (uniform_particles is preset TRUE)
    511        IF ( uniform_particles )  THEN
    512           DO  ip = nxl, nxr
    513              DO  jp = nys, nyn
    514                 DO  kp = nzb+1, nzt
    515 
    516                    n = prt_count(kp,jp,ip)
    517                    IF ( MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize  ) ==     &
    518                         MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize  )  .AND. &
    519                         MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) ==           &
    520                         MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) )  THEN
    521                       uniform_particles_l = .TRUE.
    522                    ELSE
    523                       uniform_particles_l = .FALSE.
    524                    ENDIF
    525 
    526                 ENDDO
    527              ENDDO
    528           ENDDO
    529 
    530 #if defined( __parallel )
    531           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    532           CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1,       &
    533                               MPI_LOGICAL, MPI_LAND, comm2d, ierr )
    534 #else
    535           uniform_particles = uniform_particles_l
    536 #endif
    537 
    538        ENDIF
    539 
    540 !
    541 !--    Particles will probably become none-uniform, if their size and color
    542 !--    will be determined by flow variables
    543        IF ( particle_color /= 'none'  .OR.  particle_dvrpsize /= 'none' )  THEN
    544           uniform_particles = .FALSE.
    545        ENDIF
    546 
    547 ! !kk    Not implemented aft individual particle array fort every gridcell
    548 ! !
    549 ! !--    Set the beginning of the particle tails and their age
    550 !        IF ( use_particle_tails )  THEN
    551 ! !
    552 ! !--       Choose the maximum number of tails with respect to the maximum number
    553 ! !--       of particles and skip_particles_for_tail
    554 !           maximum_number_of_tails = maximum_number_of_particles / &
    555 !                                     skip_particles_for_tail
    556 !
    557 ! !
    558 ! !--       Create a minimum number of tails in case that there is no tail
    559 ! !--       initially (otherwise, index errors will occur when adressing the
    560 ! !--       arrays below)
    561 !           IF ( maximum_number_of_tails == 0 )  maximum_number_of_tails = 10
    562 !
    563 !           ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
    564 !                     maximum_number_of_tails),                                 &
    565 !                     new_tail_id(maximum_number_of_tails),                     &
    566 !                     tail_mask(maximum_number_of_tails) )
    567 !
    568 !           particle_tail_coordinates  = 0.0_wp
    569 !           minimum_tailpoint_distance = minimum_tailpoint_distance**2
    570 !           number_of_initial_tails    = number_of_tails
    571 !
    572 !           nn = 0
    573 !           DO  n = 1, number_of_particles
    574 ! !
    575 ! !--          Only for those particles marked above with a provisional tail_id
    576 ! !--          tails will be created. Particles now get their final tail_id.
    577 !              IF ( particles(n)%tail_id /= 0 )  THEN
    578 !