Changeset 2312


Ignore:
Timestamp:
Jul 14, 2017 8:26:51 PM (4 years ago)
Author:
hoffmann
Message:

various improvements of the LCM

Location:
palm/trunk/SOURCE
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/check_parameters.f90

    r2300 r2312  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! PA0349 and PA0420 removed.
     28!
     29! 2300 2017-06-29 13:31:14Z raasch
    2730! host-specific settings and checks removed
    28 ! 
     31!
    2932! 2292 2017-06-20 09:51:42Z schwenkel
    30 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
    31 ! includes two more prognostic equations for cloud drop concentration (nc) 
    32 ! and cloud water content (qc). 
    33 ! 
     33! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     34! includes two more prognostic equations for cloud drop concentration (nc)
     35! and cloud water content (qc).
     36!
    3437! 2274 2017-06-09 13:27:48Z Giersch
    35 ! Changed error messages 
    36 ! 
     38! Changed error messages
     39!
    3740! 2271 2017-06-09 12:34:55Z sward
    3841! roughness-length check altered
    3942! Error messages fixed
    40 ! 
     43!
    4144! 2270 2017-06-09 12:18:47Z maronga
    4245! Revised numbering (removed 2 timeseries)
    43 ! 
     46!
    4447! 2259 2017-06-08 09:09:11Z gronemeier
    4548! Implemented synthetic turbulence generator
     
    5053! Doxygen comment added
    5154!
    52 ! 2248 2017-06-06 13:52:54Z sward 
     55! 2248 2017-06-06 13:52:54Z sward
    5356! Error message fixed
    5457!
    5558! 2209 2017-04-19 09:34:46Z kanani
    5659! Check for plant canopy model output
    57 ! 
     60!
    5861! 2200 2017-04-11 11:37:51Z suehring
    5962! monotonic_adjustment removed
    60 ! 
     63!
    6164! 2178 2017-03-17 11:07:39Z hellstea
    6265! Index limits for perturbations are now set also in case of nested
     
    6568! 2169 2017-03-06 18:16:35Z suehring
    6669! Bugfix, move setting for topography grid convention to init_grid
    67 ! 
     70!
    6871! 2118 2017-01-17 16:38:49Z raasch
    6972! OpenACC related parts of code removed
    70 ! 
     73!
    7174! 2088 2016-12-19 16:30:25Z suehring
    7275! Bugfix in initial salinity profile
    73 ! 
     76!
    7477! 2084 2016-12-09 15:59:42Z knoop
    7578! Removed anelastic + multigrid error message
    76 ! 
     79!
    7780! 2052 2016-11-08 15:14:59Z gronemeier
    7881! Bugfix: remove setting of default value for recycling_width
    79 ! 
     82!
    8083! 2050 2016-11-08 15:00:55Z gronemeier
    8184! Implement turbulent outflow condition
    82 ! 
     85!
    8386! 2044 2016-11-02 16:44:25Z knoop
    8487! Added error code for anelastic approximation
    85 ! 
     88!
    8689! 2042 2016-11-02 13:47:31Z suehring
    8790! Additional checks for wall_heatflux, wall_humidityflux and wall_scalarflux.
    8891! Bugfix, check for constant_scalarflux.
    89 ! 
     92!
    9093! 2040 2016-10-26 16:58:09Z gronemeier
    9194! Removed check for statistic_regions > 9.
    92 ! 
     95!
    9396! 2037 2016-10-26 11:15:40Z knoop
    9497! Anelastic approximation implemented
    95 ! 
     98!
    9699! 2031 2016-10-21 15:11:58Z knoop
    97100! renamed variable rho to rho_ocean
    98 ! 
     101!
    99102! 2026 2016-10-18 10:27:02Z suehring
    100103! Bugfix, enable output of s*2
    101 ! 
     104!
    102105! 2024 2016-10-12 16:42:37Z kanani
    103106! Added missing CASE, error message and unit for ssws*,
    104107! increased number of possible output quantities to 500.
    105 ! 
     108!
    106109! 2011 2016-09-19 17:29:57Z kanani
    107110! Flag urban_surface is now defined in module control_parameters,
    108111! changed prefix for urban surface model output to "usm_",
    109 ! added flag lsf_exception (inipar-Namelist parameter) to allow explicit 
     112! added flag lsf_exception (inipar-Namelist parameter) to allow explicit
    110113! enabling of large scale forcing together with buildings on flat terrain,
    111114! introduced control parameter varnamelength for LEN of var.
    112 ! 
     115!
    113116! 2007 2016-08-24 15:47:17Z kanani
    114117! Added checks for the urban surface model,
    115118! increased counter in DO WHILE loop over data_output (for urban surface output)
    116 ! 
     119!
    117120! 2000 2016-08-20 18:09:15Z knoop
    118121! Forced header and separation lines into 80 columns
    119 ! 
     122!
    120123! 1994 2016-08-15 09:52:21Z suehring
    121124! Add missing check for cloud_physics and cloud_droplets
    122 ! 
     125!
    123126! 1992 2016-08-12 15:14:59Z suehring
    124127! New checks for top_scalarflux
    125128! Bugfixes concerning data output of profiles in case of passive_scalar
    126 ! 
     129!
    127130! 1984 2016-08-01 14:48:05Z suehring
    128131! Bugfix: checking for bottom and top boundary condition for humidity and scalars
    129 ! 
     132!
    130133! 1972 2016-07-26 07:52:02Z maronga
    131134! Removed check of lai* (done in land surface module)
    132 ! 
     135!
    133136! 1970 2016-07-18 14:27:12Z suehring
    134137! Bugfix, check for ibc_q_b and constant_waterflux in case of land-surface scheme
    135 ! 
     138!
    136139! 1962 2016-07-13 09:16:44Z suehring
    137140! typo removed
    138 ! 
     141!
    139142! 1960 2016-07-12 16:34:24Z suehring
    140143! Separate checks and calculations for humidity and passive scalar. Therefore,
    141144! a subroutine to calculate vertical gradients is introduced, in order  to reduce
    142 ! redundance. 
     145! redundance.
    143146! Additional check - large-scale forcing is not implemented for passive scalars
    144147! so far.
    145 ! 
     148!
    146149! 1931 2016-06-10 12:06:59Z suehring
    147150! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
     
    158161! 1841 2016-04-07 19:14:06Z raasch
    159162! redundant location message removed
    160 ! 
     163!
    161164! 1833 2016-04-07 14:23:03Z raasch
    162165! check of spectra quantities moved to spectra_mod
    163 ! 
     166!
    164167! 1829 2016-04-07 12:16:29Z maronga
    165168! Bugfix: output of user defined data required reset of the string 'unit'
    166 ! 
     169!
    167170! 1826 2016-04-07 12:01:39Z maronga
    168171! Moved checks for radiation model to the respective module.
     
    172175! 1824 2016-04-07 09:55:29Z gronemeier
    173176! Check if roughness_length < dz/2. Moved location_message(finished) to the end.
    174 ! 
     177!
    175178! 1822 2016-04-07 07:49:42Z hoffmann
    176179! PALM collision kernel deleted. Collision algorithms are checked for correct
     
    186189! 1817 2016-04-06 15:44:20Z maronga
    187190! Moved checks for land_surface model to the respective module
    188 ! 
     191!
    189192! 1806 2016-04-05 18:55:35Z gronemeier
    190193! Check for recycling_yshift
     
    200203! surface presribed in the land surface model. Added output of z0q.
    201204! Syntax layout improved.
    202 ! 
     205!
    203206! 1786 2016-03-08 05:49:27Z raasch
    204207! cpp-direktives for spectra removed, check of spectra level removed
     
    217220! 1745 2016-02-05 13:06:51Z gronemeier
    218221! Bugfix: check data output intervals to be /= 0.0 in case of parallel NetCDF4
    219 ! 
     222!
    220223! 1701 2015-11-02 07:43:04Z maronga
    221224! Bugfix: definition of rad_net timeseries was missing
    222 ! 
     225!
    223226! 1691 2015-10-26 16:17:44Z maronga
    224 ! Added output of Obukhov length (ol) and radiative heating rates for RRTMG. 
     227! Added output of Obukhov length (ol) and radiative heating rates for RRTMG.
    225228! Added checks for use of radiation / lsm with topography.
    226 ! 
     229!
    227230! 1682 2015-10-07 23:56:08Z knoop
    228 ! Code annotations made doxygen readable 
    229 ! 
     231! Code annotations made doxygen readable
     232!
    230233! 1606 2015-06-29 10:43:37Z maronga
    231234! Added check for use of RRTMG without netCDF.
    232 ! 
     235!
    233236! 1587 2015-05-04 14:19:01Z maronga
    234237! Added check for set of albedos when using RRTMG
    235 ! 
     238!
    236239! 1585 2015-04-30 07:05:52Z maronga
    237240! Added support for RRTMG
     
    248251! 1555 2015-03-04 17:44:27Z maronga
    249252! Added output of r_a and r_s. Renumbering of LSM PA-messages.
    250 ! 
     253!
    251254! 1553 2015-03-03 17:33:54Z maronga
    252255! Removed check for missing soil temperature profile as default values were added.
    253 ! 
     256!
    254257! 1551 2015-03-03 14:18:16Z maronga
    255258! Added various checks for land surface and radiation model. In the course of this
    256259! action, the length of the variable var has to be increased
    257 ! 
     260!
    258261! 1504 2014-12-04 09:23:49Z maronga
    259262! Bugfix: check for large_scale forcing got lost.
    260 ! 
     263!
    261264! 1500 2014-12-03 17:42:41Z maronga
    262265! Boundary conditions changed to dirichlet for land surface model
     
    264267! 1496 2014-12-02 17:25:50Z maronga
    265268! Added checks for the land surface model
    266 ! 
     269!
    267270! 1484 2014-10-21 10:53:05Z kanani
    268271! Changes due to new module structure of the plant canopy model:
    269272!   module plant_canopy_model_mod added,
    270 !   checks regarding usage of new method for leaf area density profile 
     273!   checks regarding usage of new method for leaf area density profile
    271274!   construction added,
    272 !   lad-profile construction moved to new subroutine init_plant_canopy within 
     275!   lad-profile construction moved to new subroutine init_plant_canopy within
    273276!   the module plant_canopy_model_mod,
    274277!   drag_coefficient renamed to canopy_drag_coeff.
    275278! Missing KIND-attribute for REAL constant added
    276 ! 
     279!
    277280! 1455 2014-08-29 10:47:47Z heinze
    278 ! empty time records in volume, cross-section and masked data output prevented 
     281! empty time records in volume, cross-section and masked data output prevented
    279282! in case of non-parallel netcdf-output in restart runs
    280 ! 
     283!
    281284! 1429 2014-07-15 12:53:45Z knoop
    282285! run_description_header exended to provide ensemble_member_nr if specified
    283 ! 
     286!
    284287! 1425 2014-07-05 10:57:53Z knoop
    285288! bugfix: perturbation domain modified for parallel random number generator
    286 ! 
     289!
    287290! 1402 2014-05-09 14:25:13Z raasch
    288291! location messages modified
    289 ! 
     292!
    290293! 1400 2014-05-09 14:03:54Z knoop
    291294! Check random generator extended by option random-parallel
    292 ! 
     295!
    293296! 1384 2014-05-02 14:31:06Z raasch
    294297! location messages added
    295 ! 
     298!
    296299! 1365 2014-04-22 15:03:56Z boeske
    297300! Usage of large scale forcing for pt and q enabled:
    298 ! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q), 
    299 ! large scale subsidence (td_sub_lpt, td_sub_q) 
     301! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
     302! large scale subsidence (td_sub_lpt, td_sub_q)
    300303! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
    301304! check if use_subsidence_tendencies is used correctly
    302 ! 
     305!
    303306! 1361 2014-04-16 15:17:48Z hoffmann
    304307! PA0363 removed
    305 ! PA0362 changed 
     308! PA0362 changed
    306309!
    307310! 1359 2014-04-11 17:15:14Z hoffmann
     
    310313!
    311314! PA0084 not necessary for new particle structure
    312 ! 
     315!
    313316! 1353 2014-04-08 15:21:23Z heinze
    314317! REAL constants provided with KIND-attribute
    315 ! 
     318!
    316319! 1330 2014-03-24 17:29:32Z suehring
    317 ! In case of SGS-particle velocity advection of TKE is also allowed with 
     320! In case of SGS-particle velocity advection of TKE is also allowed with
    318321! dissipative 5th-order scheme.
    319 ! 
     322!
    320323! 1327 2014-03-21 11:00:16Z raasch
    321324! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
     
    348351! 1241 2013-10-30 11:36:58Z heinze
    349352! output for profiles of ug and vg added
    350 ! set surface_heatflux and surface_waterflux also in dependence on 
     353! set surface_heatflux and surface_waterflux also in dependence on
    351354! large_scale_forcing
    352355! checks for nudging and large scale forcing from external file
     
    395398!
    396399! 1065 2012-11-22 17:42:36Z hoffmann
    397 ! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without 
     400! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
    398401!         precipitation in order to save computational resources.
    399402!
     
    405408! - check cloud physics scheme (Kessler or Seifert and Beheng)
    406409! - plant_canopy is not allowed
    407 ! - currently, only cache loop_optimization is allowed 
     410! - currently, only cache loop_optimization is allowed
    408411! - initial profiles of nr, qr
    409412! - boundary condition of nr, qr
     
    483486!------------------------------------------------------------------------------!
    484487 SUBROUTINE check_parameters
    485  
     488
    486489
    487490    USE arrays_3d
     
    507510    USE plant_canopy_model_mod,                                                &
    508511        ONLY:  pcm_check_data_output, pcm_check_parameters, plant_canopy
    509        
     512
    510513    USE pmc_interface,                                                         &
    511514        ONLY:  cpl_id, nested_run
     
    531534    IMPLICIT NONE
    532535
    533     CHARACTER (LEN=1)   ::  sq                       !< 
    534     CHARACTER (LEN=varnamelength)  ::  var           !< 
     536    CHARACTER (LEN=1)   ::  sq                       !<
     537    CHARACTER (LEN=varnamelength)  ::  var           !<
    535538    CHARACTER (LEN=7)   ::  unit                     !<
    536     CHARACTER (LEN=8)   ::  date                     !< 
    537     CHARACTER (LEN=10)  ::  time                     !< 
     539    CHARACTER (LEN=8)   ::  date                     !<
     540    CHARACTER (LEN=10)  ::  time                     !<
    538541    CHARACTER (LEN=20)  ::  ensemble_string          !<
    539542    CHARACTER (LEN=15)  ::  nest_string              !<
    540     CHARACTER (LEN=40)  ::  coupling_string          !< 
    541     CHARACTER (LEN=100) ::  action                   !< 
    542 
    543     INTEGER(iwp) ::  i                               !< 
     543    CHARACTER (LEN=40)  ::  coupling_string          !<
     544    CHARACTER (LEN=100) ::  action                   !<
     545
     546    INTEGER(iwp) ::  i                               !<
    544547    INTEGER(iwp) ::  ilen                            !<
    545     INTEGER(iwp) ::  j                               !< 
    546     INTEGER(iwp) ::  k                               !< 
    547     INTEGER(iwp) ::  kk                              !< 
    548     INTEGER(iwp) ::  netcdf_data_format_save         !< 
     548    INTEGER(iwp) ::  j                               !<
     549    INTEGER(iwp) ::  k                               !<
     550    INTEGER(iwp) ::  kk                              !<
     551    INTEGER(iwp) ::  netcdf_data_format_save         !<
    549552    INTEGER(iwp) ::  position                        !<
    550    
    551     LOGICAL     ::  found                            !< 
    552    
    553     REAL(wp)    ::  dum                              !< 
    554     REAL(wp)    ::  gradient                         !< 
    555     REAL(wp)    ::  remote = 0.0_wp                  !< 
    556     REAL(wp)    ::  simulation_time_since_reference  !< 
     553
     554    LOGICAL     ::  found                            !<
     555
     556    REAL(wp)    ::  dum                              !<
     557    REAL(wp)    ::  gradient                         !<
     558    REAL(wp)    ::  remote = 0.0_wp                  !<
     559    REAL(wp)    ::  simulation_time_since_reference  !<
    557560
    558561
     
    568571!-- Check the coupling mode
    569572!> @todo Check if any queries for other coupling modes (e.g. precursor_ocean) are missing
    570     IF ( coupling_mode /= 'uncoupled'            .AND.  & 
    571          coupling_mode /= 'atmosphere_to_ocean'  .AND.  & 
    572          coupling_mode /= 'ocean_to_atmosphere' )  THEN     
     573    IF ( coupling_mode /= 'uncoupled'            .AND.  &
     574         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
     575         coupling_mode /= 'ocean_to_atmosphere' )  THEN
    573576       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
    574577       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
     
    590593!--    NOTE: coupled runs have not been implemented in the check_namelist_files
    591594!--    program.
    592 !--    check_namelist_files will need the following information of the other 
     595!--    check_namelist_files will need the following information of the other
    593596!--    model (atmosphere/ocean).
    594597!       dt_coupling = remote
     
    607610       ENDIF
    608611       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
    609  
     612
    610613       IF ( dt_coupling /= remote )  THEN
    611614          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
     
    620623             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
    621624                            status, ierr )
    622           ENDIF   
     625          ENDIF
    623626          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
    624      
     627
    625628          dt_coupling = MAX( dt_max, remote )
    626629          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
     
    635638          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
    636639                         status, ierr )
    637        ENDIF 
     640       ENDIF
    638641       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
    639  
     642
    640643       IF ( restart_time /= remote )  THEN
    641644          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
     
    650653          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
    651654                         status, ierr )
    652        ENDIF   
     655       ENDIF
    653656       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
    654    
     657
    655658       IF ( dt_restart /= remote )  THEN
    656659          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
     
    666669                         14, comm_inter, ierr )
    667670          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
    668                          status, ierr )   
    669        ENDIF 
     671                         status, ierr )
     672       ENDIF
    670673       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
    671    
     674
    672675       IF ( simulation_time_since_reference /= remote )  THEN
    673676          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
     
    730733             WRITE( message_string, * ) 'coupling mode "',                     &
    731734                   TRIM( coupling_mode ),                                      &
    732              '": nx+1 in ocean is not divisible by nx+1 in', & 
     735             '": nx+1 in ocean is not divisible by nx+1 in', &
    733736             ' atmosphere without remainder'
    734737             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
     
    738741             WRITE( message_string, * ) 'coupling mode "',                     &
    739742                   TRIM( coupling_mode ),                                      &
    740              '": ny+1 in ocean is not divisible by ny+1 in', & 
     743             '": ny+1 in ocean is not divisible by ny+1 in', &
    741744             ' atmosphere without remainder'
    742745
     
    763766    ENDIF
    764767    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
    765    
     768
    766769#endif
    767770
     
    932935       flux_output_mode = 'dynamic'
    933936    ENDIF
    934    
     937
    935938!
    936939!-- set the flux output units according to flux_output_mode
     
    10241027       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
    10251028    ENDIF
    1026    
    1027     IF( momentum_advec == 'ws-scheme' .AND.                                    & 
     1029
     1030    IF( momentum_advec == 'ws-scheme' .AND.                                    &
    10281031        .NOT. call_psolver_at_all_substeps  ) THEN
    10291032        message_string = 'psolver must be called at each RK3 substep when "'// &
     
    10741077    ENDIF
    10751078
    1076     IF ( use_sgs_for_particles  .AND.  curvature_solution_effects )  THEN
    1077        message_string = 'use_sgs_for_particles = .TRUE. not allowed with ' //  &
    1078                         'curvature_solution_effects = .TRUE.'
    1079        CALL message( 'check_parameters', 'PA0349', 1, 2, 0, 6, 0 )
    1080     ENDIF
    1081 
    10821079!
    10831080!-- Set LOGICAL switches to enhance performance
     
    11141111    ENDIF
    11151112!
    1116 !-- Check for proper settings for microphysics 
     1113!-- Check for proper settings for microphysics
    11171114    IF ( cloud_physics  .AND.  cloud_droplets )  THEN
    11181115       message_string = 'cloud_physics = .TRUE. is not allowed with ' //  &
     
    11411138    END SELECT
    11421139    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
    1143 
    1144 !
    1145 !-- Collision algorithms:
    1146     SELECT CASE ( TRIM( collision_algorithm ) )
    1147 
    1148        CASE ( 'all_or_nothing' )
    1149           all_or_nothing = .TRUE.
    1150 
    1151        CASE ( 'average_impact' )
    1152           average_impact = .TRUE.
    1153 
    1154        CASE DEFAULT
    1155           message_string = 'unknown collision algorithm: collision_algorithm = "' // &
    1156                            TRIM( collision_algorithm ) // '"'
    1157           CALL message( 'check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
    1158 
    1159        END SELECT
    11601140
    11611141    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     
    11861166                        ' is not allowed with conserve_volume_flow = .T.'
    11871167       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
    1188     ENDIF       
     1168    ENDIF
    11891169
    11901170
     
    12291209!-- When plant canopy model is used, peform addtional checks
    12301210    IF ( plant_canopy )  CALL pcm_check_parameters
    1231    
     1211
    12321212!
    12331213!-- Additional checks for spectra
     
    12551235                        'loop_optimization = "cache" or "vector"'
    12561236       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
    1257     ENDIF 
     1237    ENDIF
    12581238
    12591239!
     
    12881268                   i = i + 1
    12891269                ENDIF
    1290              ENDIF       
     1270             ENDIF
    12911271             IF ( gradient /= 0.0_wp )  THEN
    12921272                IF ( k /= 1 )  THEN
     
    13311311       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
    13321312          ug_vertical_gradient_level(1) = 0.0_wp
    1333        ENDIF 
     1313       ENDIF
    13341314
    13351315!
     
    13451325          vg(0) = vg_surface
    13461326          DO  k = 1, nzt+1
    1347              IF ( i < 11 )  THEN 
     1327             IF ( i < 11 )  THEN
    13481328                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
    13491329                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
     
    14821462       ENDIF
    14831463
    1484          
     1464
    14851465    ENDIF
    14861466
     
    15221502          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
    15231503        ENDIF
    1524     ENDIF   
     1504    ENDIF
    15251505
    15261506!
     
    16161596       IF ( use_ug_for_galilei_tr                    .AND.                     &
    16171597            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
    1618             ug_vertical_gradient(1) == 0.0_wp        .AND.                     & 
     1598            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
    16191599            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
    16201600            vg_vertical_gradient(1) == 0.0_wp )  THEN
     
    19131893
    19141894    ENDIF
    1915    
     1895
    19161896    IF ( passive_scalar )  THEN
    19171897
     
    19341914!
    19351915!--    A fixed scalar concentration at the top implies Dirichlet boundary
    1936 !--    condition for scalar. Hence, in this case specification of a constant 
     1916!--    condition for scalar. Hence, in this case specification of a constant
    19371917!--    scalar flux is forbidden.
    19381918       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
     
    26522632             dopr_unit(i)  = TRIM ( heatflux_output_unit )
    26532633             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
    2654              
     2634
    26552635          CASE ( 'v*pt*' )
    26562636             dopr_index(i) = 62
     
    29762956             ENDIF
    29772957
    2978  
     2958
    29792959
    29802960          CASE DEFAULT
     
    29872967                                                     unit, dopr_unit(i) )
    29882968             ENDIF
    2989              
     2969
    29902970             IF ( unit == 'illegal' )  THEN
    29912971                unit = ''
     
    32513231                                 'res passive_scalar = .TRUE.'
    32523232                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
    3253              ENDIF             
     3233             ENDIF
    32543234
    32553235             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
     
    32593239             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
    32603240             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
    3261              IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'     
     3241             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
    32623242             IF ( TRIM( var ) == 't*'     )  unit = 'K'
    32633243             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
    32643244             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
    3265              IF ( TRIM( var ) == 'z0h*'   )  unit = 'm' 
    3266              
     3245             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
     3246
    32673247          CASE DEFAULT
    32683248
    32693249             CALL lsm_check_data_output ( var, unit, i, ilen, k )
    3270  
     3250
    32713251             IF ( unit == 'illegal' )  THEN
    32723252                CALL radiation_check_data_output( var, unit, i, ilen, k )
     
    32843264                 CALL pcm_check_data_output( var, unit )
    32853265             ENDIF
    3286              
     3266
    32873267             IF ( unit == 'illegal' )  THEN
    32883268                unit = ''
     
    34083388       ENDIF
    34093389    ENDDO
    3410    
     3390
    34113391    IF ( masks < 0  .OR.  masks > max_masks )  THEN
    34123392       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
     
    34263406!
    34273407!--    Generate masks for masked data output
    3428 !--    Parallel netcdf output is not tested so far for masked data, hence 
     3408!--    Parallel netcdf output is not tested so far for masked data, hence
    34293409!--    netcdf_data_format is switched back to non-paralell output.
    34303410       netcdf_data_format_save = netcdf_data_format
     
    34363416                           '6 (parallel netCDF 4 Classic model) '//            &
    34373417                           '&are currently not supported (not yet tested) ' // &
    3438                            'for masked data.&Using respective non-parallel' // & 
     3418                           'for masked data.&Using respective non-parallel' // &
    34393419                           ' output for masked data.'
    34403420          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
     
    36293609!-- near the inflow and the perturbation area is further limited to ...(1)
    36303610!-- after the initial phase of the flow.
    3631    
     3611
    36323612    IF ( bc_lr /= 'cyclic' )  THEN
    36333613       IF ( inflow_disturbance_begin == -1 )  THEN
     
    36663646    ENDIF
    36673647
    3668     IF ( random_generator == 'random-parallel' )  THEN 
     3648    IF ( random_generator == 'random-parallel' )  THEN
    36693649       dist_nxl = nxl;  dist_nxr = nxr
    36703650       dist_nys = nys;  dist_nyn = nyn
     
    39213901
    39223902    IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
    3923        message_string = 'The usage of large scale forcing from external &'//   & 
     3903       message_string = 'The usage of large scale forcing from external &'//   &
    39243904                        'file LSF_DATA requires humidity = .T..'
    39253905       CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
     
    39273907
    39283908    IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
    3929        message_string = 'The usage of large scale forcing from external &'//   & 
     3909       message_string = 'The usage of large scale forcing from external &'//   &
    39303910                        'file LSF_DATA is not implemented for passive scalars'
    39313911       CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
     
    39343914    IF ( large_scale_forcing  .AND.  topography /= 'flat'                      &
    39353915                              .AND.  .NOT.  lsf_exception )  THEN
    3936        message_string = 'The usage of large scale forcing from external &'//   & 
     3916       message_string = 'The usage of large scale forcing from external &'//   &
    39373917                        'file LSF_DATA is not implemented for non-flat topography'
    39383918       CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
     
    39403920
    39413921    IF ( large_scale_forcing  .AND.  ocean )  THEN
    3942        message_string = 'The usage of large scale forcing from external &'//   & 
     3922       message_string = 'The usage of large scale forcing from external &'//   &
    39433923                        'file LSF_DATA is not implemented for ocean runs'
    39443924       CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
     
    39553935          do2d_yz_time_count = 0
    39563936          domask_time_count  = 0
    3957        ENDIF 
     3937       ENDIF
    39583938    ENDIF
    39593939
     
    39893969!> Check the length of data output intervals. In case of parallel NetCDF output
    39903970!> the time levels of the output files need to be fixed. Therefore setting the
    3991 !> output interval to 0.0s (usually used to output each timestep) is not 
     3971!> output interval to 0.0s (usually used to output each timestep) is not
    39923972!> possible as long as a non-fixed timestep is used.
    39933973!------------------------------------------------------------------------------!
     
    40183998
    40193999    END SUBROUTINE check_dt_do
    4020    
    4021    
    4022    
    4023    
     4000
     4001
     4002
     4003
    40244004!------------------------------------------------------------------------------!
    40254005! Description:
     
    40344014
    40354015
    4036        IMPLICIT NONE   
    4037    
     4016       IMPLICIT NONE
     4017
    40384018       INTEGER(iwp) ::  i     !< counter
    40394019       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
    4040        
    4041        REAL(wp)     ::  bc_t_val      !< model top gradient       
    4042        REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity 
     4020
     4021       REAL(wp)     ::  bc_t_val      !< model top gradient
     4022       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
    40434023       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
    40444024
    4045        
     4025
    40464026       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
    40474027       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
    40484028       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
    4049        
     4029
    40504030       i = 1
    40514031       gradient = 0.0_wp
     
    41094089
    41104090       ENDIF
    4111        
     4091
    41124092!
    41134093!--    In case of no given gradients, choose zero gradient conditions
     
    41184098!--    Store gradient at the top boundary for possible Neumann boundary condition
    41194099       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
    4120    
     4100
    41214101    END SUBROUTINE init_vertical_profiles
    4122    
    4123    
    4124      
     4102
     4103
     4104
    41254105!------------------------------------------------------------------------------!
    41264106! Description:
     
    41294109!------------------------------------------------------------------------------!
    41304110
    4131     SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t ) 
    4132 
    4133 
    4134        IMPLICIT NONE   
    4135    
    4136        CHARACTER (LEN=1)   ::  sq                       !< 
     4111    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
     4112
     4113
     4114       IMPLICIT NONE
     4115
     4116       CHARACTER (LEN=1)   ::  sq                       !<
    41374117       CHARACTER (LEN=*)   ::  bc_b
    41384118       CHARACTER (LEN=*)   ::  bc_t
     
    41714151          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
    41724152       ENDIF
    4173        
    4174    
    4175     END SUBROUTINE set_bc_scalars   
     4153
     4154
     4155    END SUBROUTINE set_bc_scalars
    41764156
    41774157
     
    41804160! Description:
    41814161! ------------
    4182 !> Check for consistent settings of bottom boundary conditions for humidity 
     4162!> Check for consistent settings of bottom boundary conditions for humidity
    41834163!> and scalars.
    41844164!------------------------------------------------------------------------------!
     
    41894169
    41904170
    4191        IMPLICIT NONE   
    4192    
    4193        CHARACTER (LEN=1)   ::  sq                       !< 
     4171       IMPLICIT NONE
     4172
     4173       CHARACTER (LEN=1)   ::  sq                       !<
    41944174       CHARACTER (LEN=*)   ::  bc_b
    41954175       CHARACTER (LEN=*)   ::  err_nr_1
    41964176       CHARACTER (LEN=*)   ::  err_nr_2
    4197        
     4177
    41984178       INTEGER(iwp)        ::  ibc_b
    4199        
     4179
    42004180       LOGICAL             ::  constant_flux
    4201        
     4181
    42024182       REAL(wp)            ::  surface_initial_change
    4203  
     4183
    42044184!
    42054185!--    A given surface value implies Dirichlet boundary condition for
     
    42204200                 surface_initial_change
    42214201          CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
    4222        ENDIF 
    4223        
    4224    
    4225     END SUBROUTINE check_bc_scalars   
    4226    
    4227    
     4202       ENDIF
     4203
     4204
     4205    END SUBROUTINE check_bc_scalars
     4206
     4207
    42284208
    42294209 END SUBROUTINE check_parameters
  • palm/trunk/SOURCE/data_output_ptseries.f90

    r2101 r2312  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! SGS velocities also possible for curvature_solution_effects = .TRUE.
    2728!
    2829! 2000 2016-08-20 18:09:15Z knoop
     
    3637!
    3738! 1682 2015-10-07 23:56:08Z knoop
    38 ! Code annotations made doxygen readable 
    39 ! 
     39! Code annotations made doxygen readable
     40!
    4041! 1359 2014-04-11 17:15:14Z hoffmann
    41 ! New particle structure integrated. 
    42 ! 
     42! New particle structure integrated.
     43!
    4344! 1353 2014-04-08 15:21:23Z heinze
    4445! REAL constants provided with KIND-attribute
    45 ! 
     46!
    4647! 1327 2014-03-21 11:00:16Z raasch
    4748! -netcdf output queries
     
    7576!------------------------------------------------------------------------------!
    7677 SUBROUTINE data_output_ptseries
    77  
     78
    7879
    7980    USE control_parameters,                                                    &
     
    9798
    9899    USE particle_attributes,                                                   &
    99         ONLY:  curvature_solution_effects, grid_particles, number_of_particles,&
    100                number_of_particle_groups, particles, prt_count
     100        ONLY:  grid_particles, number_of_particles, number_of_particle_groups, &
     101               particles, prt_count
    101102
    102103    USE pegrid
     
    164165                   pts_value_l(0,7)  = pts_value_l(0,7) + particles(n)%speed_y  ! mean v
    165166                   pts_value_l(0,8)  = pts_value_l(0,8) + particles(n)%speed_z  ! mean w
    166                    IF ( .NOT. curvature_solution_effects )  THEN
    167                       pts_value_l(0,9)  = pts_value_l(0,9)  + particles(n)%rvar1 ! mean sgsu
    168                       pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
    169                       pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
    170                    ENDIF
     167                   pts_value_l(0,9)  = pts_value_l(0,9)  + particles(n)%rvar1 ! mean sgsu
     168                   pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
     169                   pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
    171170                   IF ( particles(n)%speed_z > 0.0_wp )  THEN
    172171                      pts_value_l(0,12) = pts_value_l(0,12) + 1.0_wp  ! # of upward moving prts
     
    198197                      pts_value_l(jg,7)  = pts_value_l(jg,7) + particles(n)%speed_y
    199198                      pts_value_l(jg,8)  = pts_value_l(jg,8) + particles(n)%speed_z
    200                       IF ( .NOT. curvature_solution_effects )  THEN
    201                          pts_value_l(jg,9)  = pts_value_l(jg,9)  + particles(n)%rvar1
    202                          pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
    203                          pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
    204                       ENDIF
     199                      pts_value_l(jg,9)  = pts_value_l(jg,9)  + particles(n)%rvar1
     200                      pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
     201                      pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
    205202                      IF ( particles(n)%speed_z > 0.0_wp )  THEN
    206203                         pts_value_l(jg,12) = pts_value_l(jg,12) + 1.0_wp
     
    299296                pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
    300297                                                          pts_value(0,8) )**2   ! w*2
    301                 IF ( .NOT. curvature_solution_effects )  THEN
    302                    pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
    303                                                              pts_value(0,9) )**2   ! u"2
    304                    pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
    305                                                              pts_value(0,10) )**2  ! v"2
    306                    pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
    307                                                              pts_value(0,11) )**2  ! w"2
    308                 ENDIF
     298                pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
     299                                                          pts_value(0,9) )**2   ! u"2
     300                pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
     301                                                          pts_value(0,10) )**2  ! v"2
     302                pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
     303                                                          pts_value(0,11) )**2  ! w"2
    309304!
    310305!--             Repeat the same for the respective particle group
     
    324319                   pts_value_l(jg,25) = pts_value_l(jg,25) + ( particles(n)%speed_z - &
    325320                                                             pts_value(jg,8) )**2
    326                    IF ( .NOT. curvature_solution_effects )  THEN
    327                       pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
    328                                                                 pts_value(jg,9) )**2
    329                       pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
    330                                                                 pts_value(jg,10) )**2
    331                       pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
    332                                                                 pts_value(jg,11) )**2
    333                    ENDIF
     321                   pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
     322                                                             pts_value(jg,9) )**2
     323                   pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
     324                                                             pts_value(jg,10) )**2
     325                   pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
     326                                                             pts_value(jg,11) )**2
    334327                ENDIF
    335328
  • palm/trunk/SOURCE/lpm_droplet_collision.f90

    r2274 r2312  
    2020! Current revisions:
    2121! ------------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Consideration of aerosol mass during collision. Average impact algorithm has
     28! been removed.
     29!
     30! 2274 2017-06-09 13:27:48Z Giersch
    2731! Changed error messages
    2832!
     
    3943!
    4044! 1860 2016-04-13 13:21:28Z hoffmann
    41 ! Interpolation of dissipation rate adjusted to more reasonable values. 
     45! Interpolation of dissipation rate adjusted to more reasonable values.
    4246!
    4347! 1822 2016-04-07 07:49:42Z hoffmann
     
    4953!
    5054! 1682 2015-10-07 23:56:08Z knoop
    51 ! Code annotations made doxygen readable 
    52 ! 
     55! Code annotations made doxygen readable
     56!
    5357! 1359 2014-04-11 17:15:14Z hoffmann
    54 ! New particle structure integrated. 
     58! New particle structure integrated.
    5559! Kind definition added to all floating point numbers.
    5660!
     
    96100!>              the modification of the relative velocity between the droplets,
    97101!>              the effect of preferential concentration, and the enhancement of
    98 !>              collision efficiencies. 
     102!>              collision efficiencies.
    99103!------------------------------------------------------------------------------!
    100104 SUBROUTINE lpm_droplet_collision (i,j,k)
    101  
    102 
    103105
    104106    USE arrays_3d,                                                             &
     
    126128
    127129    USE particle_attributes,                                                   &
    128         ONLY:  all_or_nothing, average_impact, dissipation_classes,            &
    129                hall_kernel, iran_part, number_of_particles, particles,         &
    130                particle_type, prt_count, use_kernel_tables, wang_kernel
     130        ONLY:  curvature_solution_effects, dissipation_classes, hall_kernel,   &
     131               iran_part, number_of_particles, particles, particle_type,       &
     132               prt_count, rho_s, use_kernel_tables, wang_kernel
    131133
    132134    USE random_function_mod,                                                   &
     
    150152    REAL(wp) ::  epsilon                 !< dissipation rate
    151153    REAL(wp) ::  factor_volume_to_mass   !< 4.0 / 3.0 * pi * rho_l
    152     REAL(wp) ::  xm                      !< mean mass of droplet m
    153     REAL(wp) ::  xn                      !< mean mass of droplet n
    154 
    155     REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight  !< weighting factor
    156     REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass    !< total mass of super droplet
     154    REAL(wp) ::  xm                      !< droplet mass of super-droplet m
     155    REAL(wp) ::  xn                      !< droplet mass of super-droplet n
     156    REAL(wp) ::  xsm                     !< aerosol mass of super-droplet m
     157    REAL(wp) ::  xsn                     !< aerosol mass of super-droplet n
     158
     159    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight    !< weighting factor
     160    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass      !< total mass of super droplet
     161    REAL(wp), DIMENSION(:), ALLOCATABLE ::  aero_mass !< total aerosol mass of super droplet
    157162
    158163    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
    159164
    160165    number_of_particles   = prt_count(k,j,i)
    161     factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l 
    162     ddV                   = 1 / ( dx * dy * dz )
     166    factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l
     167    ddV                   = 1.0_wp / ( dx * dy * dz )
    163168!
    164169!-- Collision requires at least one super droplet inside the box
    165170    IF ( number_of_particles > 0 )  THEN
    166171
    167 !
    168 !--    Now apply the different kernels
    169172       IF ( use_kernel_tables )  THEN
    170173!
    171174!--       Fast method with pre-calculated collection kernels for
    172175!--       discrete radius- and dissipation-classes.
    173 !--
    174 !--       Determine dissipation class index of this gridbox
    175176          IF ( wang_kernel )  THEN
    176177             eclass = INT( diss(k,j,i) * 1.0E4_wp / 600.0_wp * &
     
    180181             epsilon = 0.0_wp
    181182          ENDIF
     183
    182184          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
    183185             eclass = 0   ! Hall kernel is used
     
    186188          ENDIF
    187189
    188 !
    189 !--       Droplet collision are calculated using collision-coalescence
    190 !--       formulation proposed by Wang (see PALM documentation)
    191 !--       Temporary fields for total mass of super-droplet and weighting factors
    192 !--       are allocated.
    193           ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
    194 
    195           mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
    196                                           particles(1:number_of_particles)%radius**3     * &
    197                                           factor_volume_to_mass
    198           weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
    199 
    200           IF ( average_impact )  THEN  ! select collision algorithm
    201 
    202              DO  n = 1, number_of_particles
    203 
    204                 rclass_l = particles(n)%class
    205                 xn       = mass(n) / weight(n)
    206 
    207                 DO  m = n, number_of_particles
    208 
    209                    rclass_s = particles(m)%class
    210                    xm       = mass(m) / weight(m)
    211 
    212                    IF ( xm .LT. xn )  THEN
    213                      
    214 !
    215 !--                   Particle n collects smaller particle m
    216                       collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
    217                                                weight(n) * ddV * dt_3d
    218 
    219                       mass(n)   = mass(n)   + mass(m)   * collection_probability
    220                       weight(m) = weight(m) - weight(m) * collection_probability
    221                       mass(m)   = mass(m)   - mass(m)   * collection_probability
    222                    ELSEIF ( xm .GT. xn )  THEN
    223 !
    224 !--                   Particle m collects smaller particle n
    225                       collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
    226                                                weight(m) * ddV * dt_3d
    227 
    228                       mass(m)   = mass(m)   + mass(n)   * collection_probability
    229                       weight(n) = weight(n) - weight(n) * collection_probability
    230                       mass(n)   = mass(n)   - mass(n)   * collection_probability
    231                    ELSE
    232 !
    233 !--                   Same-size collections. If n = m, weight is reduced by the
    234 !--                   number of possible same-size collections; the total mass
    235 !--                   is not changed during same-size collection.
    236 !--                   Same-size collections of different
    237 !--                   particles ( n /= m ) are treated as same-size collections
    238 !--                   of ONE partilce with weight = weight(n) + weight(m) and
    239 !--                   mass = mass(n) + mass(m).
    240 !--                   Accordingly, each particle loses the same number of
    241 !--                   droplets to the other particle, but this has no effect on
    242 !--                   total mass mass, since the exchanged droplets have the
    243 !--                   same radius.
    244 
    245 !--                   Note: For m = n this equation is an approximation only
    246 !--                   valid for weight >> 1 (which is usually the case). The
    247 !--                   approximation is weight(n)-1 = weight(n).
    248                       weight(n) = weight(n) - 0.5_wp * weight(n) *                &
    249                                               ckernel(rclass_l,rclass_s,eclass) * &
    250                                               weight(m) * ddV * dt_3d
    251                       IF ( n .NE. m )  THEN
    252                          weight(m) = weight(m) - 0.5_wp * weight(m) *                &
    253                                                  ckernel(rclass_l,rclass_s,eclass) * &
    254                                                  weight(n) * ddV * dt_3d
    255                       ENDIF
    256                    ENDIF
    257 
    258                 ENDDO
    259 
    260                 ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
    261 
    262              ENDDO
    263 
    264           ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
    265 
    266              DO  n = 1, number_of_particles
    267 
    268                 rclass_l = particles(n)%class
    269                 xn       = mass(n) / weight(n) ! mean mass of droplet n
    270 
    271                 DO  m = n, number_of_particles
    272 
    273                    rclass_s = particles(m)%class
    274                    xm = mass(m) / weight(m) ! mean mass of droplet m
    275 
    276                    IF ( weight(n) .LT. weight(m) )  THEN
    277 !
    278 !--                   Particle n collects weight(n) droplets of particle m 
    279                       collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
    280                                                weight(m) * ddV * dt_3d
    281 
    282                       IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    283                          mass(n)   = mass(n)   + weight(n) * xm
    284                          weight(m) = weight(m) - weight(n)
    285                          mass(m)   = mass(m)   - weight(n) * xm
    286                       ENDIF
    287 
    288                    ELSEIF ( weight(m) .LT. weight(n) )  THEN
    289 !
    290 !--                   Particle m collects weight(m) droplets of particle n 
    291                       collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
    292                                                weight(n) * ddV * dt_3d
    293 
    294                       IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    295                          mass(m)   = mass(m)   + weight(m) * xn
    296                          weight(n) = weight(n) - weight(m)
    297                          mass(n)   = mass(n)   - weight(m) * xn
    298                       ENDIF
    299                    ELSE
    300 !
    301 !--                   Collisions of particles of the same weighting factor.
    302 !--                   Particle n collects 1/2 weight(n) droplets of particle m,
    303 !--                   particle m collects 1/2 weight(m) droplets of particle n.
    304 !--                   The total mass mass changes accordingly.
    305 !--                   If n = m, the first half of the droplets coalesces with the
    306 !--                   second half of the droplets; mass is unchanged because
    307 !--                   xm = xn for n = m.
    308 
    309 !--                   Note: For m = n this equation is an approximation only
    310 !--                   valid for weight >> 1 (which is usually the case). The
    311 !--                   approximation is weight(n)-1 = weight(n).
    312                       collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
    313                                                weight(n) * ddV * dt_3d
    314 
    315                       IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    316                          mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
    317                          mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
    318                          weight(n) = weight(n) - 0.5_wp * weight(m)
    319                          weight(m) = weight(n)
    320                       ENDIF
    321                    ENDIF
    322 
    323                 ENDDO
    324 
    325                 ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
    326 
    327              ENDDO
    328 
    329           ENDIF
    330 
    331 
    332 
    333 
    334           IF ( ANY(weight < 0.0_wp) )  THEN
    335                 WRITE( message_string, * ) 'negative weighting factor'
    336                 CALL message( 'lpm_droplet_collision', 'PA0028',      &
    337                                2, 2, -1, 6, 1 )
    338           ENDIF
    339 
    340           particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
    341                                                       ( weight(1:number_of_particles) &
    342                                                         * factor_volume_to_mass       &
    343                                                       )                               &
    344                                                     )**0.33333333333333_wp
    345 
    346           particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
    347 
    348           DEALLOCATE(weight, mass)
    349 
    350        ELSEIF ( .NOT. use_kernel_tables )  THEN
    351 !
    352 !--       Collection kernels are calculated for every new
     190       ELSE
     191!
     192!--       Collection kernels are re-calculated for every new
    353193!--       grid box. First, allocate memory for kernel table.
    354194!--       Third dimension is 1, because table is re-calculated for
     
    359199!--       the kernel is based on the previous time step
    360200          CALL recalculate_kernel( i, j, k )
    361 !
    362 !--       Droplet collision are calculated using collision-coalescence
    363 !--       formulation proposed by Wang (see PALM documentation)
    364 !--       Temporary fields for total mass of super-droplet and weighting factors
    365 !--       are allocated.
    366           ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
    367 
    368           mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
    369                                         particles(1:number_of_particles)%radius**3     * &
    370                                         factor_volume_to_mass
    371 
    372           weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
    373 
    374           IF ( average_impact )  THEN  ! select collision algorithm
    375 
    376              DO  n = 1, number_of_particles
    377 
    378                 xn = mass(n) / weight(n) ! mean mass of droplet n
    379 
    380                 DO  m = n, number_of_particles
    381 
    382                    xm = mass(m) / weight(m) !mean mass of droplet m
    383 
    384                    IF ( xm .LT. xn )  THEN
    385 !
    386 !--                   Particle n collects smaller particle m
    387                       collection_probability = ckernel(n,m,1) * weight(n) *    &
    388                                                ddV * dt_3d
    389 
    390                       mass(n)   = mass(n)   + mass(m)   * collection_probability
    391                       weight(m) = weight(m) - weight(m) * collection_probability
    392                       mass(m)   = mass(m)   - mass(m)   * collection_probability
    393                    ELSEIF ( xm .GT. xn )  THEN
    394 !
    395 !--                   Particle m collects smaller particle n
    396                       collection_probability = ckernel(n,m,1) * weight(m) *    &
    397                                                ddV * dt_3d
    398 
    399                       mass(m)   = mass(m)   + mass(n)   * collection_probability
    400                       weight(n) = weight(n) - weight(n) * collection_probability
    401                       mass(n)   = mass(n)   - mass(n)   * collection_probability
    402                    ELSE
    403 !
    404 !--                   Same-size collections. If n = m, weight is reduced by the
    405 !--                   number of possible same-size collections; the total mass
    406 !--                   mass is not changed during same-size collection.
    407 !--                   Same-size collections of different
    408 !--                   particles ( n /= m ) are treated as same-size collections
    409 !--                   of ONE partilce with weight = weight(n) + weight(m) and
    410 !--                   mass = mass(n) + mass(m).
    411 !--                   Accordingly, each particle loses the same number of
    412 !--                   droplets to the other particle, but this has no effect on
    413 !--                   total mass mass, since the exchanged droplets have the
    414 !--                   same radius.
     201
     202       ENDIF
     203!
     204!--    Temporary fields for total mass of super-droplet, aerosol mass, and
     205!--    weighting factor are allocated.
     206       ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
     207       IF ( curvature_solution_effects )  ALLOCATE(aero_mass(1:number_of_particles))
     208
     209       mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
     210                                       particles(1:number_of_particles)%radius**3     * &
     211                                       factor_volume_to_mass
     212
     213       weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
     214
     215       IF ( curvature_solution_effects )  THEN
     216          aero_mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
     217                                             particles(1:number_of_particles)%aux1**3       * &
     218                                             4.0 / 3.0 * pi * rho_s
     219       ENDIF
     220!
     221!--    Calculate collision/coalescence
     222       DO  n = 1, number_of_particles
     223
     224          DO  m = n, number_of_particles
     225!
     226!--          For collisions, the weighting factor of at least one super-droplet
     227!--          needs to be larger or equal to one.
     228             IF ( MIN( weight(n), weight(m) ) .LT. 1.0 )  CYCLE
     229!
     230!--          Get mass of individual droplets (aerosols)
     231             xn = mass(n) / weight(n)
     232             xm = mass(m) / weight(m)
     233             IF ( curvature_solution_effects )  THEN
     234                xsn = aero_mass(n) / weight(n)
     235                xsm = aero_mass(m) / weight(m)
     236             ENDIF
     237!
     238!--          Probability that the necessary collisions take place
     239             IF ( use_kernel_tables )  THEN
     240                rclass_l = particles(n)%class
     241                rclass_s = particles(m)%class
     242
     243                collection_probability  = MAX( weight(n), weight(m) ) *     &
     244                                          ckernel(rclass_l,rclass_s,eclass) * ddV * dt_3d
     245             ELSE
     246                collection_probability  = MAX( weight(n), weight(m) ) *     &
     247                                          ckernel(n,m,1) * ddV * dt_3d
     248             ENDIF
     249!
     250!--          Calculate the number of collections and consider multiple collections.
     251!--          (Accordingly, p_crit will be 0.0, 1.0, 2.0, ...)
     252             IF ( collection_probability - FLOOR(collection_probability)    &
     253                  .GT. random_function( iran_part ) )  THEN
     254                collection_probability = FLOOR(collection_probability) + 1.0_wp
     255             ELSE
     256                collection_probability = FLOOR(collection_probability)
     257             ENDIF
     258
     259             IF ( collection_probability .GT. 0.0 )  THEN
     260!
     261!--             Super-droplet n collects droplets of super-droplet m
     262                IF ( weight(n) .LT. weight(m) )  THEN
     263
     264                   mass(n)   = mass(n)   + weight(n) * xm * collection_probability
     265                   weight(m) = weight(m) - weight(n)      * collection_probability
     266                   mass(m)   = mass(m)   - weight(n) * xm * collection_probability
     267                   IF ( curvature_solution_effects )  THEN
     268                      aero_mass(n) = aero_mass(n) + weight(n) * xsm * collection_probability
     269                      aero_mass(m) = aero_mass(m) - weight(n) * xsm * collection_probability
     270                   ENDIF
     271
     272                ELSEIF ( weight(m) .LT. weight(n) )  THEN
     273
     274                   mass(m)   = mass(m)   + weight(m) * xn * collection_probability
     275                   weight(n) = weight(n) - weight(m)      * collection_probability
     276                   mass(n)   = mass(n)   - weight(m) * xn * collection_probability
     277                   IF ( curvature_solution_effects )  THEN
     278                      aero_mass(m) = aero_mass(m) + weight(m) * xsn * collection_probability
     279                      aero_mass(n) = aero_mass(n) - weight(m) * xsn * collection_probability
     280                   ENDIF
     281
     282                ELSE
     283!
     284!--                Collisions of particles of the same weighting factor.
     285!--                Particle n collects 1/2 weight(n) droplets of particle m,
     286!--                particle m collects 1/2 weight(m) droplets of particle n.
     287!--                The total mass mass changes accordingly.
     288!--                If n = m, the first half of the droplets coalesces with the
     289!--                second half of the droplets; mass is unchanged because
     290!--                xm = xn for n = m.
    415291!--
    416 !--                   Note: For m = n this equation is an approximation only
    417 !--                   valid for weight >> 1 (which is usually the case). The
    418 !--                   approximation is weight(n)-1 = weight(n).
    419                       weight(n) = weight(n) - 0.5_wp * weight(n) *             &
    420                                               ckernel(n,m,1) * weight(m) *     &
    421                                               ddV * dt_3d
    422                       IF ( n .NE. m )  THEN
    423                          weight(m) = weight(m) - 0.5_wp * weight(m) *          &
    424                                                  ckernel(n,m,1) * weight(n) *  &
    425                                                  ddV * dt_3d
    426                       ENDIF
     292!--                Note: For m = n this equation is an approximation only
     293!--                valid for weight >> 1 (which is usually the case). The
     294!--                approximation is weight(n)-1 = weight(n).
     295                   mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
     296                   mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
     297                   IF ( curvature_solution_effects )  THEN
     298                      aero_mass(n) = aero_mass(n) + 0.5_wp * weight(n) * ( xsm - xsn )
     299                      aero_mass(m) = aero_mass(m) + 0.5_wp * weight(m) * ( xsn - xsm )
    427300                   ENDIF
    428 
    429 
    430                 ENDDO
    431 
    432                 ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
    433 
    434              ENDDO
    435 
    436           ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
    437 
    438              DO  n = 1, number_of_particles
    439 
    440                 xn = mass(n) / weight(n) ! mean mass of droplet n
    441 
    442                 DO  m = n, number_of_particles
    443 
    444                    xm = mass(m) / weight(m) !mean mass of droplet m
    445 
    446                    IF ( weight(n) .LT. weight(m) )  THEN
    447 !
    448 !--                   Particle n collects smaller particle m
    449                       collection_probability = ckernel(n,m,1) * weight(m) *    &
    450                                                ddV * dt_3d
    451 
    452                       IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    453                          mass(n)   = mass(n)   + weight(n) * xm 
    454                          weight(m) = weight(m) - weight(n)
    455                          mass(m)   = mass(m)   - weight(n) * xm
    456                       ENDIF
    457 
    458                    ELSEIF ( weight(m) .LT. weight(n) )  THEN
    459 !
    460 !--                   Particle m collects smaller particle n
    461                       collection_probability = ckernel(n,m,1) * weight(n) *    &
    462                                                ddV * dt_3d
    463 
    464                       IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    465                          mass(m)   = mass(m)   + weight(m) * xn
    466                          weight(n) = weight(n) - weight(m)
    467                          mass(n)   = mass(n)   - weight(m) * xn
    468                       ENDIF
    469                    ELSE
    470 !
    471 !--                   Collisions of particles of the same weighting factor.
    472 !--                   Particle n collects 1/2 weight(n) droplets of particle m,
    473 !--                   particle m collects 1/2 weight(m) droplets of particle n.
    474 !--                   The total mass mass changes accordingly.
    475 !--                   If n = m, the first half of the droplets coalesces with the
    476 !--                   second half of the droplets; mass is unchanged because
    477 !--                   xm = xn for n = m.
    478 !--
    479 !--                   Note: For m = n this equation is an approximation only
    480 !--                   valid for weight >> 1 (which is usually the case). The
    481 !--                   approximation is weight(n)-1 = weight(n).
    482                       collection_probability = ckernel(n,m,1) * weight(n) *    &
    483                                                ddV * dt_3d
    484 
    485                       IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    486                          mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
    487                          mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
    488                          weight(n) = weight(n) - 0.5_wp * weight(m)
    489                          weight(m) = weight(n)
    490                       ENDIF
    491                    ENDIF
    492 
    493 
    494                 ENDDO
    495 
    496                 ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
    497 
    498              ENDDO
    499 
    500           ENDIF
    501 
    502           IF ( ANY(weight < 0.0_wp) )  THEN
    503                 WRITE( message_string, * ) 'negative weighting factor'
    504                 CALL message( 'lpm_droplet_collision', 'PA0028',      &
    505                                2, 2, -1, 6, 1 )
    506           ENDIF
    507 
    508           particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
    509                                                       ( weight(1:number_of_particles) &
    510                                                         * factor_volume_to_mass       &
    511                                                       )                               &
    512                                                     )**0.33333333333333_wp
    513 
    514           particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
    515 
    516           DEALLOCATE( weight, mass, ckernel )
    517 
    518        ENDIF
    519  
     301                   weight(n) = weight(n) - 0.5_wp * weight(m)
     302                   weight(m) = weight(n)
     303
     304                ENDIF
     305
     306             ENDIF
     307
     308          ENDDO
     309
     310          ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     311
     312       ENDDO
     313
     314       IF ( ANY(weight < 0.0_wp) )  THEN
     315             WRITE( message_string, * ) 'negative weighting factor'
     316             CALL message( 'lpm_droplet_collision', 'PA0028',      &
     317                            2, 2, -1, 6, 1 )
     318       ENDIF
     319
     320       particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
     321                                                   ( weight(1:number_of_particles) &
     322                                                     * factor_volume_to_mass       &
     323                                                   )                               &
     324                                                 )**0.33333333333333_wp
     325
     326       IF ( curvature_solution_effects )  THEN
     327          particles(1:number_of_particles)%aux1 = ( aero_mass(1:number_of_particles) / &
     328                                                    ( weight(1:number_of_particles)    &
     329                                                      * 4.0_wp / 3.0_wp * pi * rho_s   &
     330                                                    )                                  &
     331                                                  )**0.33333333333333_wp
     332       ENDIF
     333
     334       particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
     335
     336       DEALLOCATE( weight, mass )
     337       IF ( curvature_solution_effects )  DEALLOCATE( aero_mass )
     338       IF ( .NOT. use_kernel_tables )  DEALLOCATE( ckernel )
     339
    520340!
    521341!--    Check if LWC is conserved during collision process
     
    532352
    533353    ENDIF
    534  
     354
    535355    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
    536356
  • palm/trunk/SOURCE/lpm_droplet_condensation.f90

    r2101 r2312  
    2020! Current revisions:
    2121! ------------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Rosenbrock scheme improved. Gas-kinetic effect added.
    2728!
    2829! 2000 2016-08-20 18:09:15Z knoop
    2930! Forced header and separation lines into 80 columns
    30 ! 
     31!
    3132! 1890 2016-04-22 08:52:11Z hoffmann
    3233! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more
    3334! than 40 iterations to find a sufficient time setp, the model is not aborted.
    34 ! This might lead to small erros. But they can be assumend as negligible, since 
    35 ! the maximum timesetp of the Rosenbrock method after 40 iterations will be 
    36 ! smaller than 10^-16 s. 
    37 ! 
     35! This might lead to small erros. But they can be assumend as negligible, since
     36! the maximum timesetp of the Rosenbrock method after 40 iterations will be
     37! smaller than 10^-16 s.
     38!
    3839! 1871 2016-04-15 11:46:09Z hoffmann
    3940! Initialization of aerosols added.
     
    5354!
    5455! 1682 2015-10-07 23:56:08Z knoop
    55 ! Code annotations made doxygen readable 
    56 ! 
     56! Code annotations made doxygen readable
     57!
    5758! 1359 2014-04-11 17:15:14Z hoffmann
    58 ! New particle structure integrated. 
     59! New particle structure integrated.
    5960! Kind definition added to all floating point numbers.
    6061!
    6162! 1346 2014-03-27 13:18:20Z heinze
    62 ! Bugfix: REAL constants provided with KIND-attribute especially in call of 
     63! Bugfix: REAL constants provided with KIND-attribute especially in call of
    6364! intrinsic function like MAX, MIN, SIGN
    6465!
     
    108109!------------------------------------------------------------------------------!
    109110 SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
    110  
     111
    111112
    112113    USE arrays_3d,                                                             &
    113         ONLY:  hyp, pt, q,  ql_c, ql_v
     114        ONLY:  hyp, pt, q, ql_c, ql_v
    114115
    115116    USE cloud_parameters,                                                      &
     
    139140               use_kernel_tables, vanthoff, wang_kernel
    140141
    141 
    142142    IMPLICIT NONE
    143143
    144144    INTEGER(iwp) :: ip                         !<
    145     INTEGER(iwp) :: internal_timestep_count    !<
    146145    INTEGER(iwp) :: jp                         !<
    147     INTEGER(iwp) :: jtry                       !<
    148146    INTEGER(iwp) :: kp                         !<
    149147    INTEGER(iwp) :: n                          !<
    150     INTEGER(iwp) :: ros_count                  !<
    151  
    152     INTEGER(iwp), PARAMETER ::  maxtry = 40    !<
    153 
    154     LOGICAL ::  repeat                         !<
    155 
    156     REAL(wp) ::  aa                            !<
     148
    157149    REAL(wp) ::  afactor                       !< curvature effects
    158150    REAL(wp) ::  arg                           !<
     
    161153    REAL(wp) ::  delta_r                       !<
    162154    REAL(wp) ::  diameter                      !< diameter of cloud droplets
    163     REAL(wp) ::  diff_coeff_v                  !< diffusivity for water vapor
     155    REAL(wp) ::  diff_coeff                    !< diffusivity for water vapor
    164156    REAL(wp) ::  drdt                          !<
    165     REAL(wp) ::  drdt_ini                      !<
    166157    REAL(wp) ::  dt_ros                        !<
    167     REAL(wp) ::  dt_ros_next                   !<
    168158    REAL(wp) ::  dt_ros_sum                    !<
    169     REAL(wp) ::  dt_ros_sum_ini                !<
    170159    REAL(wp) ::  d2rdtdr                       !<
    171     REAL(wp) ::  errmax                        !<
    172160    REAL(wp) ::  e_a                           !< current vapor pressure
    173161    REAL(wp) ::  e_s                           !< current saturation vapor pressure
    174     REAL(wp) ::  err_ros                       !<
    175     REAL(wp) ::  g1                            !<
    176     REAL(wp) ::  g2                            !<
    177     REAL(wp) ::  g3                            !<
    178     REAL(wp) ::  g4                            !<
    179     REAL(wp) ::  r_ros                         !<
    180     REAL(wp) ::  r_ros_ini                     !<
    181     REAL(wp) ::  sigma                         !<
    182     REAL(wp) ::  thermal_conductivity_v        !< thermal conductivity for water
     162    REAL(wp) ::  error                         !< local truncation error in Rosenbrock
     163    REAL(wp) ::  k1                            !<
     164    REAL(wp) ::  k2                            !<
     165    REAL(wp) ::  r_err                         !< First order estimate of Rosenbrock radius
     166    REAL(wp) ::  r_ros                         !< Rosenbrock radius
     167    REAL(wp) ::  r_ros_ini                     !< initial Rosenbrock radius
     168    REAL(wp) ::  r0                            !< gas-kinetic lengthscale
     169    REAL(wp) ::  sigma                         !< surface tension of water
     170    REAL(wp) ::  thermal_conductivity          !< thermal conductivity for water
    183171    REAL(wp) ::  t_int                         !< temperature
    184172    REAL(wp) ::  w_s                           !< terminal velocity of droplets
    185     REAL(wp) ::  re_p                          !<
    186 
    187 !
    188 !-- Parameters for Rosenbrock method
    189     REAL(wp), PARAMETER ::  a21 = 2.0_wp               !<
    190     REAL(wp), PARAMETER ::  a31 = 48.0_wp / 25.0_wp    !<
    191     REAL(wp), PARAMETER ::  a32 = 6.0_wp / 25.0_wp     !<
    192     REAL(wp), PARAMETER ::  b1 = 19.0_wp / 9.0_wp      !<
    193     REAL(wp), PARAMETER ::  b2 = 0.5_wp                !<
    194     REAL(wp), PARAMETER ::  b3 = 25.0_wp / 108.0_wp    !<
    195     REAL(wp), PARAMETER ::  b4 = 125.0_wp / 108.0_wp   !<
    196     REAL(wp), PARAMETER ::  c21 = -8.0_wp              !<
    197     REAL(wp), PARAMETER ::  c31 = 372.0_wp / 25.0_wp   !<
    198     REAL(wp), PARAMETER ::  c32 = 12.0_wp / 5.0_wp     !<
    199     REAL(wp), PARAMETER ::  c41 = -112.0_wp / 125.0_wp !<
    200     REAL(wp), PARAMETER ::  c42 = -54.0_wp / 125.0_wp  !<
    201     REAL(wp), PARAMETER ::  c43 = -2.0_wp / 5.0_wp     !<
    202     REAL(wp), PARAMETER ::  errcon = 0.1296_wp         !<
    203     REAL(wp), PARAMETER ::  e1 = 17.0_wp / 54.0_wp     !<
    204     REAL(wp), PARAMETER ::  e2 = 7.0_wp / 36.0_wp      !<
    205     REAL(wp), PARAMETER ::  e3 = 0.0_wp                !<
    206     REAL(wp), PARAMETER ::  e4 = 125.0_wp / 108.0_wp   !<
    207     REAL(wp), PARAMETER ::  eps_ros = 1.0E-3_wp        !< accuracy of Rosenbrock method
    208     REAL(wp), PARAMETER ::  gam = 0.5_wp               !<
    209     REAL(wp), PARAMETER ::  grow = 1.5_wp              !<
    210     REAL(wp), PARAMETER ::  pgrow = -0.25_wp           !<
    211     REAL(wp), PARAMETER ::  pshrnk = -1.0_wp /3.0_wp   !<
    212     REAL(wp), PARAMETER ::  shrnk = 0.5_wp             !<
    213 
     173    REAL(wp) ::  re_p                          !< particle Reynolds number
     174!
     175!-- Parameters for Rosenbrock method (see Verwer et al., 1999)
     176    REAL(wp), PARAMETER :: prec = 1.0E-3_wp     !< precision of Rosenbrock solution
     177    REAL(wp), PARAMETER :: q_increase = 1.5_wp  !< increase factor in timestep
     178    REAL(wp), PARAMETER :: q_decrease = 0.9_wp  !< decrease factor in timestep
     179    REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
    214180!
    215181!-- Parameters for terminal velocity
     
    224190    REAL(wp), DIMENSION(number_of_particles) ::  new_r                  !<
    225191
    226 
    227 
    228192    CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
    229193
    230194!
    231 !-- Calculate temperature, saturation vapor pressure and current vapor pressure
     195!-- Absolute temperature
    232196    t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
    233     e_s   = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) )
    234     e_a   = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp )
    235 !
    236 !-- Thermal conductivity for water (from Rogers and Yau, Table 7.1),
    237 !-- diffusivity for water vapor (after Hall und Pruppacher, 1976)
    238     thermal_conductivity_v = 7.94048E-05_wp * t_int + 0.00227011_wp
    239     diff_coeff_v           = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * &
    240                              ( 101325.0_wp / hyp(kp) )
    241 !
    242 !-- Calculate effects of heat conductivity and diffusion of water vapor on the
    243 !-- condensation/evaporation process (typically known as 1.0 / (F_k + F_d) )
    244     ddenom  = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff_v ) +        &
     197!
     198!-- Saturation vapor pressure (Eq. 10 in Bolton, 1980)
     199    e_s = 611.2_wp * EXP( 17.62_wp * ( t_int - 273.15_wp ) / ( t_int - 29.65_wp ) )
     200!
     201!-- Current vapor pressure
     202    e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + 0.622_wp )
     203!
     204!-- Thermal conductivity for water (from Rogers and Yau, Table 7.1)
     205    thermal_conductivity = 7.94048E-05_wp * t_int + 0.00227011_wp
     206!
     207!-- Moldecular diffusivity of water vapor in air (Hall und Pruppacher, 1976)
     208    diff_coeff           = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * &
     209                           ( 101325.0_wp / hyp(kp) )
     210!
     211!-- Lengthscale for gas-kinetic effects (from Mordy, 1959, p. 23):
     212    r0 = diff_coeff / 0.036_wp * SQRT( 2.0_wp * pi / ( r_v * t_int ) )
     213!
     214!-- Calculate effects of heat conductivity and diffusion of water vapor on the
     215!-- diffusional growth process (usually known as 1.0 / (F_k + F_d) )
     216    ddenom  = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff ) +          &
    245217                         ( l_v / ( r_v * t_int ) - 1.0_wp ) * rho_l *          &
    246                          l_v / ( thermal_conductivity_v * t_int )              &
     218                         l_v / ( thermal_conductivity * t_int )                &
    247219                       )
    248 
    249220    new_r = 0.0_wp
    250 
    251221!
    252222!-- Determine ventilation effect on evaporation of large drops
     
    264234          ENDIF
    265235!
    266 !--       First calculate droplet's Reynolds number
     236!--       Calculate droplet's Reynolds number
    267237          re_p = 2.0_wp * particles(n)%radius * w_s / molecular_viscosity
    268238!
     
    275245       ELSE
    276246!
    277 !--       For small droplets or in supersaturated environments, the ventilation 
     247!--       For small droplets or in supersaturated environments, the ventilation
    278248!--       effect does not play a role
    279249          ventilation_effect(n) = 1.0_wp
     
    281251    ENDDO
    282252
    283 !
    284 !-- Use analytic model for condensational growth
    285253    IF( .NOT. curvature_solution_effects ) then
     254!
     255!--    Use analytic model for diffusional growth including gas-kinetic
     256!--    effects (Mordy, 1959) but without the impact of aerosols.
    286257       DO  n = 1, number_of_particles
    287           arg = particles(n)%radius**2 + 2.0_wp * dt_3d * ddenom *            &
    288                                          ventilation_effect(n) *               &
    289                                          ( e_a / e_s - 1.0_wp )
    290           arg = MAX( arg, 1.0E-16_wp )
    291           new_r(n) = SQRT( arg )
     258          arg      = ( particles(n)%radius + r0 )**2 + 2.0_wp * dt_3d * ddenom * &
     259                                                       ventilation_effect(n) *   &
     260                                                       ( e_a / e_s - 1.0_wp )
     261          arg      = MAX( arg, ( 0.01E-6 + r0 )**2 )
     262          new_r(n) = SQRT( arg ) - r0
    292263       ENDDO
    293     ENDIF
    294 
    295 !
    296 !-- If selected, use numerical solution of the condensational growth
    297 !-- equation (e.g., for studying the activation of aerosols).
    298 !-- Curvature and solutions effects are included in growth equation.
    299 !-- Change in Radius is calculated with a 4th-order Rosenbrock method
    300 !-- for stiff o.d.e's with monitoring local truncation error to adjust
    301 !-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
    302     DO  n = 1, number_of_particles
    303        IF ( curvature_solution_effects )  THEN
    304 
    305           ros_count = 0
    306           repeat = .TRUE.
    307 !
    308 !--       Carry out the Rosenbrock algorithm. In case of unreasonable results
    309 !--       the switch "repeat" will be set true and the algorithm will be carried
    310 !--       out again with the internal time step set to its initial (small) value.
    311 !--       Unreasonable results may occur if the external conditions, especially
    312 !--       the supersaturation, has significantly changed compared to the last
    313 !--       PALM timestep.
    314           DO WHILE ( repeat )
    315 
    316              repeat = .FALSE.
    317 !
    318 !--          Curvature effect (afactor) with surface tension parameterization
    319 !--          by Straka (2009)
    320              sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
    321              afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
    322 !
    323 !--          Solute effect (bfactor)
    324              bfactor = vanthoff * rho_s * particles(n)%rvar2**3 *              &
    325                        molecular_weight_of_water /                             &
    326                        ( rho_l * molecular_weight_of_solute )
    327 
    328              r_ros = particles(n)%radius
    329              dt_ros_sum  = 0.0_wp      ! internal integrated time (s)
    330              internal_timestep_count = 0
    331 !
    332 !--          Take internal time step values from the end of last PALM time step
    333              dt_ros_next = particles(n)%rvar1
    334 
    335 !
    336 !--          Internal time step should not be > 1.0E-2 and < 1.0E-6
    337              IF ( dt_ros_next > 1.0E-2_wp )  THEN
    338                 dt_ros_next = 1.0E-2_wp
    339              ELSEIF ( dt_ros_next < 1.0E-6_wp )  THEN
    340                 dt_ros_next = 1.0E-6_wp
    341              ENDIF
    342 
    343 !
    344 !--          If calculation of Rosenbrock method is repeated due to unreasonalble
    345 !--          results during previous try the initial internal time step has to be
    346 !--          reduced
    347              IF ( ros_count > 1 )  THEN
    348                 dt_ros_next = dt_ros_next * 0.1_wp
    349              ELSEIF ( ros_count > 5 )  THEN
    350 !
    351 !--             Prevent creation of infinite loop
    352                 message_string = 'ros_count > 5 in Rosenbrock method'
    353                 CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, &
    354                                0, 6, 0 )
    355              ENDIF
    356 
    357 !
    358 !--          Internal time step must not be larger than PALM time step
    359              dt_ros = MIN( dt_ros_next, dt_3d )
    360 
    361 !
    362 !--          Integrate growth equation in time unless PALM time step is reached
    363              DO WHILE ( dt_ros_sum < dt_3d )
    364 
    365                 internal_timestep_count = internal_timestep_count + 1
    366 
    367 !
    368 !--             Derivative at starting value
    369                 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
     264
     265    ELSE
     266!
     267!--    Integrate the diffusional growth including gas-kinetic (Mordy, 1959),
     268!--    as well as curvature and solute effects (e.g., Köhler, 1936).
     269!
     270!--    Curvature effect (afactor) with surface tension (sigma) by Straka (2009)
     271       sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
     272!
     273!--    Solute effect (afactor)
     274       afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
     275
     276       DO  n = 1, number_of_particles
     277!
     278!--       Solute effect (bfactor)
     279          bfactor = vanthoff * rho_s * particles(n)%aux1**3 *                    &
     280                    molecular_weight_of_water / ( rho_l * molecular_weight_of_solute )
     281
     282          dt_ros     = particles(n)%aux2  ! use previously stored Rosenbrock timestep
     283          dt_ros_sum = 0.0_wp
     284
     285          r_ros     = particles(n)%radius  ! initialize Rosenbrock particle radius
     286          r_ros_ini = r_ros
     287!
     288!--       Integrate growth equation using a 2nd-order Rosenbrock method
     289!--       (see Verwer et al., 1999, Eq. (3.2)). The Rosenbrock method adjusts
     290!--       its with internal timestep to minimize the local truncation error.
     291          DO WHILE ( dt_ros_sum < dt_3d )
     292
     293             dt_ros = MIN( dt_ros, dt_3d - dt_ros_sum )
     294
     295             DO
     296
     297                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -    &
    370298                                                          afactor / r_ros +    &
    371299                                                          bfactor / r_ros**3   &
    372                                                         ) / r_ros
    373 
    374                 drdt_ini       = drdt
    375                 dt_ros_sum_ini = dt_ros_sum
    376                 r_ros_ini      = r_ros
    377 
    378 !
    379 !--             Calculate radial derivative of dr/dt
    380                 d2rdtdr = ddenom * ventilation_effect(n) *                     &
    381                                        ( ( 1.0_wp - e_a / e_s ) / r_ros**2 +   &
    382                                          2.0_wp * afactor / r_ros**3 -         &
    383                                          4.0_wp * bfactor / r_ros**5           &
    384                                        )
    385 !
    386 !--             Adjust stepsize unless required accuracy is reached
    387                 DO  jtry = 1, maxtry+1
    388 
    389                    IF ( jtry == maxtry+1 )  THEN
    390                       message_string = 'maxtry > 40 in Rosenbrock method'
    391                       CALL message( 'lpm_droplet_condensation', 'PA0347', 0,   &
    392                                     1, 0, 6, 0 )
    393                    ENDIF
    394 
    395                    aa    = 1.0_wp / ( gam * dt_ros ) - d2rdtdr
    396                    g1    = drdt_ini / aa
    397                    r_ros = r_ros_ini + a21 * g1
    398                    drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
    399                                                               afactor / r_ros +    &
    400                                                               bfactor / r_ros**3   &
    401                                                             ) / r_ros
    402 
    403                    g2    = ( drdt + c21 * g1 / dt_ros )&
    404                            / aa
    405                    r_ros = r_ros_ini + a31 * g1 + a32 * g2
    406                    drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
    407                                                               afactor / r_ros +    &
    408                                                               bfactor / r_ros**3   &
    409                                                             ) / r_ros
    410 
    411                    g3    = ( drdt +  &
    412                              ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
    413                    g4    = ( drdt +  &
    414                              ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
    415                    r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
    416 
    417                    dt_ros_sum = dt_ros_sum_ini + dt_ros
    418 
    419                    IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
    420                       message_string = 'zero stepsize in Rosenbrock method'
    421                       CALL message( 'lpm_droplet_condensation', 'PA0348', 2,   &
    422                                     2, 0, 6, 0 )
    423                    ENDIF
    424 !
    425 !--                Calculate error
    426                    err_ros = e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4
    427                    errmax  = 0.0_wp
    428                    errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
    429 !
    430 !--                Leave loop if accuracy is sufficient, otherwise try again
    431 !--                with a reduced stepsize
    432                    IF ( errmax <= 1.0_wp )  THEN
    433                       EXIT
    434                    ELSE
    435                       dt_ros = MAX( ABS( 0.9_wp * dt_ros * errmax**pshrnk ),   &
    436                                     shrnk * ABS( dt_ros ) )
    437                    ENDIF
    438 
    439                 ENDDO  ! loop for stepsize adjustment
    440 
    441 !
    442 !--             Calculate next internal time step
    443                 IF ( errmax > errcon )  THEN
    444                    dt_ros_next = 0.9_wp * dt_ros * errmax**pgrow
     300                                                        ) / ( r_ros + r0 )
     301
     302                d2rdtdr = -ddenom * ventilation_effect(n) * (                  &
     303                                                (e_a / e_s - 1.0) * r_ros**4 - &
     304                                                afactor * r0 * r_ros**2 -      &
     305                                                2.0 * afactor * r_ros**3 +     &
     306                                                3.0 * bfactor * r0 +           &
     307                                                4.0 * bfactor * r_ros          &
     308                                                            )                  &
     309                          / ( r_ros**4 * ( r_ros + r0 )**2 )
     310
     311                k1    = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )
     312
     313                r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1)
     314                r_err = r_ros
     315
     316                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -   &
     317                                                           afactor / r_ros +   &
     318                                                           bfactor / r_ros**3  &
     319                                                         ) / ( r_ros + r0 )
     320
     321                k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / &
     322                     ( 1.0 - dt_ros * gamma * d2rdtdr )
     323
     324                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1)
     325   !
     326   !--          Check error of the solution, and reduce dt_ros if necessary.
     327                error = ABS(r_err - r_ros) / r_ros
     328                IF ( error .GT. prec )  THEN
     329                   dt_ros = SQRT( q_decrease * prec / error ) * dt_ros
     330                   r_ros  = r_ros_ini
    445331                ELSE
    446                    dt_ros_next = grow * dt_ros
     332                   dt_ros_sum = dt_ros_sum + dt_ros
     333                   dt_ros     = q_increase * dt_ros
     334                   r_ros_ini  = r_ros
     335                   EXIT
    447336                ENDIF
    448337
    449 !
    450 !--             Estimated time step is reduced if the PALM time step is exceeded
    451                 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
    452                    dt_ros = dt_3d - dt_ros_sum
    453                 ELSE
    454                    dt_ros = dt_ros_next
    455                 ENDIF
    456 
    457              ENDDO
    458 !
    459 !--          Store internal time step value for next PALM step
    460              particles(n)%rvar1 = dt_ros_next
    461 
    462 !
    463 !--          Radius should not fall below 1E-8 because Rosenbrock method may
    464 !--          lead to errors otherwise
    465              new_r(n) = MAX( r_ros, particles(n)%rvar2 )
    466 !
    467 !--          Check if calculated droplet radius change is reasonable since in
    468 !--          case of droplet evaporation the Rosenbrock method may lead to
    469 !--          secondary solutions which are physically unlikely.
    470 !--          Due to the solution effect the droplets may grow for relative
    471 !--          humidities below 100%, but change of radius should not be too
    472 !--          large. In case of unreasonable droplet growth the Rosenbrock
    473 !--          method is recalculated using a smaller initial time step.
    474 !--          Limiting values are tested for droplets down to 1.0E-7
    475              IF ( new_r(n) - particles(n)%radius >= 3.0E-7_wp  .AND.  &
    476                   e_a / e_s < 0.97_wp )  THEN
    477                 ros_count = ros_count + 1
    478                 repeat = .TRUE.
    479              ENDIF
    480 
    481           ENDDO    ! Rosenbrock method
    482 
    483        ENDIF
    484 
    485        delta_r = new_r(n) - particles(n)%radius
    486 
    487 !
    488 !--    Sum up the change in volume of liquid water for the respective grid
    489 !--    volume (this is needed later in lpm_calc_liquid_water_content for
    490 !--    calculating the release of latent heat)
     338             END DO
     339
     340          END DO !Rosenbrock loop
     341!
     342!--       Store new particle radius
     343          new_r(n) = r_ros
     344!
     345!--       Store internal time step value for next PALM step
     346          particles(n)%aux2 = dt_ros
     347
     348       ENDDO !Particle loop
     349
     350    ENDIF
     351
     352    DO  n = 1, number_of_particles
     353!
     354!--    Sum up the change in liquid water for the respective grid
     355!--    box for the computation of the release/depletion of water vapor
     356!--    and heat.
    491357       ql_c(kp,jp,ip) = ql_c(kp,jp,ip) + particles(n)%weight_factor *          &
    492358                                   rho_l * 1.33333333_wp * pi *                &
    493359                                   ( new_r(n)**3 - particles(n)%radius**3 ) /  &
    494360                                   ( rho_surface * dx * dy * dz )
     361!
     362!--    Check if the increase in liqid water is not too big. If this is the case,
     363!--    the model timestep might be too long.
    495364       IF ( ql_c(kp,jp,ip) > 100.0_wp )  THEN
    496365          WRITE( message_string, * ) 'k=',kp,' j=',jp,' i=',ip,      &
     
    499368          CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
    500369       ENDIF
    501 
    502 !
    503 !--    Change the droplet radius
    504        IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp  .AND.        &
    505            new_r(n) < 0.0_wp )  THEN
     370!
     371!--    Check if the change in the droplet radius is not too big. If this is the
     372!--    case, the model timestep might be too long.
     373       delta_r = new_r(n) - particles(n)%radius
     374       IF ( delta_r < 0.0_wp  .AND. new_r(n) < 0.0_wp )  THEN
    506375          WRITE( message_string, * ) '#1 k=',kp,' j=',jp,' i=',ip,    &
    507376                       ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,      &
     
    510379          CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
    511380       ENDIF
    512 
    513381!
    514382!--    Sum up the total volume of liquid water (needed below for
    515383!--    re-calculating the weighting factors)
    516384       ql_v(kp,jp,ip) = ql_v(kp,jp,ip) + particles(n)%weight_factor * new_r(n)**3
    517 
    518        particles(n)%radius = new_r(n)
    519 
    520385!
    521386!--    Determine radius class of the particle needed for collision
    522        IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables )     &
    523        THEN
     387       IF ( use_kernel_tables )  THEN
    524388          particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) /           &
    525389                               ( rclass_ubound - rclass_lbound ) *             &
     
    528392          particles(n)%class = MAX( particles(n)%class, 1 )
    529393       ENDIF
     394 !
     395 !--   Store new radius to particle features
     396       particles(n)%radius = new_r(n)
    530397
    531398    ENDDO
  • palm/trunk/SOURCE/lpm_init.f90

    r2305 r2312  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Extended particle data type. Aerosol initialization improved.
     28!
     29! 2305 2017-07-06 11:18:47Z hoffmann
    2730! Improved calculation of particle IDs.
    28 ! 
     31!
    2932! 2274 2017-06-09 13:27:48Z Giersch
    3033!  Changed error messages
    31 ! 
     34!
    3235! 2265 2017-06-08 16:58:28Z schwenkel
    3336! Unused variables removed.
    34 ! 
     37!
    3538! 2263 2017-06-08 14:59:01Z schwenkel
    3639! Implemented splitting and merging algorithm
    37 ! 
     40!
    3841! 2233 2017-05-30 18:08:54Z suehring
    3942!
    4043! 2232 2017-05-30 17:47:52Z suehring
    4144! Adjustments according to new topography realization
    42 ! 
    43 ! 
     45!
     46!
    4447! 2223 2017-05-15 16:38:09Z suehring
    4548! Add check for particle release at model top
    46 ! 
     49!
    4750! 2182 2017-03-17 14:27:40Z schwenkel
    4851! Added parameters for simplified particle initialization.
    4952!
    5053! 2122 2017-01-18 12:22:54Z hoffmann
    51 ! Improved initialization of equilibrium aerosol radii 
     54! Improved initialization of equilibrium aerosol radii
    5255! Calculation of particle ID
    5356!
    5457! 2000 2016-08-20 18:09:15Z knoop
    5558! Forced header and separation lines into 80 columns
    56 ! 
     59!
    5760! 2016-06-09 16:25:25Z suehring
    58 ! Bugfix in determining initial particle height and grid index in case of 
     61! Bugfix in determining initial particle height and grid index in case of
    5962! seed_follows_topography.
    60 ! Bugfix concerning random positions, ensure that particles do not move more 
     63! Bugfix concerning random positions, ensure that particles do not move more
    6164! than one grid length.
    6265! Bugfix logarithmic interpolation.
     
    6467!
    6568! 1890 2016-04-22 08:52:11Z hoffmann
    66 ! Initialization of aerosol equilibrium radius not possible in supersaturated 
     69! Initialization of aerosol equilibrium radius not possible in supersaturated
    6770! environments. Therefore, a maximum supersaturation of -1 % is assumed during
    6871! initialization.
     
    7073! 1873 2016-04-18 14:50:06Z maronga
    7174! Module renamed (removed _mod
    72 ! 
     75!
    7376! 1871 2016-04-15 11:46:09Z hoffmann
    7477! Initialization of aerosols added.
     
    8689! netcdf module added
    8790!
    88 ! 1725 2015-11-17 13:01:51Z hoffmann 
    89 ! Bugfix: Processor-dependent seed for random function is generated before it is 
     91! 1725 2015-11-17 13:01:51Z hoffmann
     92! Bugfix: Processor-dependent seed for random function is generated before it is
    9093! used.
    9194!
     
    97100!
    98101! 1682 2015-10-07 23:56:08Z knoop
    99 ! Code annotations made doxygen readable 
     102! Code annotations made doxygen readable
    100103!
    101104! 1575 2015-03-27 09:56:27Z raasch
     
    103106!
    104107! 1359 2014-04-11 17:15:14Z hoffmann
    105 ! New particle structure integrated. 
     108! New particle structure integrated.
    106109! Kind definition added to all floating point numbers.
    107110! lpm_init changed form a subroutine to a module.
    108 ! 
     111!
    109112! 1327 2014-03-21 11:00:16Z raasch
    110113! -netcdf_output
     
    123126!
    124127! 1314 2014-03-14 18:25:17Z suehring
    125 ! Vertical logarithmic interpolation of horizontal particle speed for particles 
     128! Vertical logarithmic interpolation of horizontal particle speed for particles
    126129! between roughness height and first vertical grid level.
    127130!
     
    146149!
    147150! 667 2010-12-23 12:06:00Z suehring/gryschka
    148 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation 
     151! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
    149152! of arrays.
    150153!
     
    159162!------------------------------------------------------------------------------!
    160163 MODULE lpm_init_mod
    161  
     164
    162165    USE, INTRINSIC ::  ISO_C_BINDING
    163166
     
    167170    USE control_parameters,                                                    &
    168171        ONLY:  cloud_droplets, constant_flux_layer, current_timestep_number,   &
    169                dz, initializing_actions, message_string, ocean, simulated_time
     172               dt_3d, dz, initializing_actions, message_string, ocean,         &
     173               simulated_time
    170174
    171175    USE grid_variables,                                                        &
     
    197201                particles, particle_advection_start, particle_groups,          &
    198202                particle_groups_type, particles_per_point,                     &
    199                 particle_type, pdx, pdy, pdz,  prt_count, psb, psl, psn, psr,  & 
    200                 pss, pst, radius, random_start_position,                       & 
     203                particle_type, pdx, pdy, pdz,  prt_count, psb, psl, psn, psr,  &
     204                pss, pst, radius, random_start_position,                       &
    201205                read_particles_from_restartfile, seed_follows_topography,      &
    202206                sgs_wf_part, sort_count, splitting_function, splitting_mode,   &
     
    247251    INTEGER(iwp) ::  k                           !<
    248252
    249     REAL(wp) ::  div                             !< 
     253    REAL(wp) ::  div                             !<
    250254    REAL(wp) ::  height_int                      !<
    251255    REAL(wp) ::  height_p                        !<
     
    286290!-- Check if downward-facing walls exist. This case, reflection boundary
    287291!-- conditions (as well as subgrid-scale velocities) may do not work
    288 !-- propably (not realized so far). 
     292!-- propably (not realized so far).
    289293    IF ( surf_def_h(1)%ns >= 1 )  THEN
    290        WRITE( message_string, * ) 'Overhanging topography do not work '// &   
    291                                   'with particles' 
     294       WRITE( message_string, * ) 'Overhanging topography do not work '// &
     295                                  'with particles'
    292296       CALL message( 'lpm_init', 'PA0212', 0, 1, 0, 6, 0 )
    293297
    294     ENDIF 
     298    ENDIF
    295299
    296300!
     
    310314!-- If number_particles_per_gridbox is set, the parametres pdx, pdy and pdz are
    311315!-- calculated diagnostically. Therfore an isotropic distribution is prescribed.
    312     IF ( number_particles_per_gridbox /= -1 .AND.   & 
     316    IF ( number_particles_per_gridbox /= -1 .AND.   &
    313317         number_particles_per_gridbox >= 1 )    THEN
    314        pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  & 
     318       pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  &
    315319             REAL(number_particles_per_gridbox))**0.3333333_wp
    316320!
    317 !--    Ensure a smooth value (two significant digits) of distance between 
     321!--    Ensure a smooth value (two significant digits) of distance between
    318322!--    particles (pdx, pdy, pdz).
    319323       div = 1000.0_wp
     
    340344
    341345!
    342 !-- Allocate arrays required for calculating particle SGS velocities. 
     346!-- Allocate arrays required for calculating particle SGS velocities.
    343347!-- Initialize prefactor required for stoachastic Weil equation.
    344348    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
     
    347351                 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    348352
    349        sgs_wf_part = 1.0_wp / 3.0_wp   
     353       sgs_wf_part = 1.0_wp / 3.0_wp
    350354    ENDIF
    351355
    352356!
    353 !-- Allocate array required for logarithmic vertical interpolation of 
     357!-- Allocate array required for logarithmic vertical interpolation of
    354358!-- horizontal particle velocities between the surface and the first vertical
    355 !-- grid level. In order to avoid repeated CPU cost-intensive CALLS of 
    356 !-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for 
    357 !-- several heights. Splitting into 20 sublayers turned out to be sufficient. 
     359!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
     360!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
     361!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
    358362!-- To obtain exact height levels of particles, linear interpolation is applied
    359363!-- (see lpm_advec.f90).
    360364    IF ( constant_flux_layer )  THEN
    361        
    362        ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 
     365
     366       ALLOCATE ( log_z_z0(0:number_of_sublayers) )
    363367       z_p = zu(nzb+1) - zw(nzb)
    364368
    365369!
    366370!--    Calculate horizontal mean value of z0 used for logartihmic
    367 !--    interpolation. Note: this is not exact for heterogeneous z0. 
    368 !--    However, sensitivity studies showed that the effect is 
    369 !--    negligible. 
     371!--    interpolation. Note: this is not exact for heterogeneous z0.
     372!--    However, sensitivity studies showed that the effect is
     373!--    negligible.
    370374       z0_av_local  = SUM( surf_def_h(0)%z0 ) + SUM( surf_lsm_h%z0 ) +         &
    371375                      SUM( surf_usm_h%z0 )
     
    401405!-- Check boundary condition and set internal variables
    402406    SELECT CASE ( bc_par_b )
    403    
     407
    404408       CASE ( 'absorb' )
    405409          ibc_par_b = 1
     
    407411       CASE ( 'reflect' )
    408412          ibc_par_b = 2
    409          
     413
    410414       CASE DEFAULT
    411415          WRITE( message_string, * )  'unknown boundary condition ',   &
    412416                                       'bc_par_b = "', TRIM( bc_par_b ), '"'
    413417          CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
    414          
     418
    415419    END SELECT
    416420    SELECT CASE ( bc_par_t )
    417    
     421
    418422       CASE ( 'absorb' )
    419423          ibc_par_t = 1
     
    421425       CASE ( 'reflect' )
    422426          ibc_par_t = 2
    423          
     427
    424428       CASE DEFAULT
    425429          WRITE( message_string, * ) 'unknown boundary condition ',   &
    426430                                     'bc_par_t = "', TRIM( bc_par_t ), '"'
    427431          CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
    428          
     432
    429433    END SELECT
    430434    SELECT CASE ( bc_par_lr )
     
    438442       CASE ( 'reflect' )
    439443          ibc_par_lr = 2
    440          
     444
    441445       CASE DEFAULT
    442446          WRITE( message_string, * ) 'unknown boundary condition ',   &
    443447                                     'bc_par_lr = "', TRIM( bc_par_lr ), '"'
    444448          CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
    445          
     449
    446450    END SELECT
    447451    SELECT CASE ( bc_par_ns )
     
    455459       CASE ( 'reflect' )
    456460          ibc_par_ns = 2
    457          
     461
    458462       CASE DEFAULT
    459463          WRITE( message_string, * ) 'unknown boundary condition ',   &
    460464                                     'bc_par_ns = "', TRIM( bc_par_ns ), '"'
    461465          CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
    462          
     466
    463467    END SELECT
    464468    SELECT CASE ( splitting_mode )
    465    
     469
    466470       CASE ( 'const' )
    467471          i_splitting_mode = 1
     
    472476       CASE ( 'gb_av' )
    473477          i_splitting_mode = 3
    474          
     478
    475479       CASE DEFAULT
    476480          WRITE( message_string, * )  'unknown splitting condition ',   &
    477481                                       'splitting_mode = "', TRIM( splitting_mode ), '"'
    478482          CALL message( 'lpm_init', 'PA0146', 1, 2, 0, 6, 0 )
    479          
     483
    480484    END SELECT
    481485    SELECT CASE ( splitting_function )
    482    
     486
    483487       CASE ( 'gamma' )
    484488          isf = 1
     
    489493       CASE ( 'exp' )
    490494          isf = 3
    491          
     495
    492496       CASE DEFAULT
    493497          WRITE( message_string, * )  'unknown splitting function ',   &
    494498                                       'splitting_function = "', TRIM( splitting_function ), '"'
    495499          CALL message( 'lpm_init', 'PA0147', 1, 2, 0, 6, 0 )
    496          
     500
    497501    END SELECT
    498    
     502
    499503
    500504!
     
    524528
    525529!
    526 !--    initialize counter for particle IDs 
     530!--    initialize counter for particle IDs
    527531       grid_particles%id_counter = 1
    528532
     
    534538                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    535539                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    536                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
     540                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    537541                                      0, 0, 0_idp, .FALSE., -1 )
    538542
     
    577581          CALL check_open( 80 )
    578582          WRITE ( 80, 8000 )  current_timestep_number, simulated_time,         &
    579                               number_of_particles                             
     583                              number_of_particles
    580584          CALL close_file( 80 )
    581585       ENDIF
     
    584588
    585589!
    586 !-- To avoid programm abort, assign particles array to the local version of 
     590!-- To avoid programm abort, assign particles array to the local version of
    587591!-- first grid cell
    588592    number_of_particles = prt_count(nzb+1,nys,nxl)
     
    608612
    609613    USE particle_attributes,                                                   &
    610         ONLY: deleted_particles, monodisperse_aerosols
     614        ONLY: deleted_particles
    611615
    612616    IMPLICIT  NONE
     
    631635    LOGICAL                    ::  first_stride !< flag for initialization
    632636
    633     REAL(wp)                   ::  pos_x      !< increment for particle position in x     
    634     REAL(wp)                   ::  pos_y      !< increment for particle position in y 
    635     REAL(wp)                   ::  pos_z      !< increment for particle position in z     
     637    REAL(wp)                   ::  pos_x      !< increment for particle position in x
     638    REAL(wp)                   ::  pos_y      !< increment for particle position in y
     639    REAL(wp)                   ::  pos_z      !< increment for particle position in z
    636640    REAL(wp)                   ::  rand_contr !< dummy argument for random position
    637641
     
    652656!--    Calculate initial_weighting_factor diagnostically
    653657       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
    654           initial_weighting_factor =  number_concentration * 1.0E6_wp *             & 
    655                                       pdx(1) * pdy(1) * pdz(1) 
     658          initial_weighting_factor =  number_concentration * 1.0E6_wp *             &
     659                                      pdx(1) * pdy(1) * pdz(1)
    656660       END IF
    657661
     
    677681               xloop: DO WHILE ( pos_x <= psr(i) )
    678682
    679                          IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            & 
     683                         IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            &
    680684                              pos_x <  ( nxr + 0.5_wp ) * dx )  THEN
    681685
     
    690694                               tmp_particle%age_m         = 0.0_wp
    691695                               tmp_particle%dt_sum        = 0.0_wp
    692                                tmp_particle%user          = 0.0_wp !unused, free for the user
    693696                               tmp_particle%e_m           = 0.0_wp
    694                                IF ( curvature_solution_effects )  THEN
    695 !
    696 !--                               Initial values (internal timesteps, derivative)
    697 !--                               for Rosenbrock method
    698                                   tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
    699                                   tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
    700                                   tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
    701                                ELSE
    702 !
    703 !--                               Initial values for SGS velocities
    704                                   tmp_particle%rvar1      = 0.0_wp
    705                                   tmp_particle%rvar2      = 0.0_wp
    706                                   tmp_particle%rvar3      = 0.0_wp
    707                                ENDIF
     697                               tmp_particle%rvar1         = 0.0_wp
     698                               tmp_particle%rvar2         = 0.0_wp
     699                               tmp_particle%rvar3         = 0.0_wp
    708700                               tmp_particle%speed_x       = 0.0_wp
    709701                               tmp_particle%speed_y       = 0.0_wp
     
    712704                               tmp_particle%origin_y      = pos_y
    713705                               tmp_particle%origin_z      = pos_z
     706                               IF ( curvature_solution_effects )  THEN
     707                                  tmp_particle%aux1      = 0.0_wp    ! dry aerosol radius
     708                                  tmp_particle%aux2      = dt_3d     ! last Rosenbrock timestep
     709                               ELSE
     710                                  tmp_particle%aux1      = 0.0_wp    ! free to use
     711                                  tmp_particle%aux2      = 0.0_wp    ! free to use
     712                               ENDIF
    714713                               tmp_particle%radius        = particle_groups(i)%radius
    715714                               tmp_particle%weight_factor = initial_weighting_factor
     
    726725!
    727726!--                            Determine surface level. Therefore, check for
    728 !--                            upward-facing wall on w-grid. MAXLOC will return 
     727!--                            upward-facing wall on w-grid. MAXLOC will return
    729728!--                            the index of the lowest upward-facing wall.
    730729                               k_surf = MAXLOC(                                &
     
    748747                                  ENDIF
    749748!
    750 !--                            Skip particle release if particle position is 
     749!--                            Skip particle release if particle position is
    751750!--                            below surface, or within topography in case
    752751!--                            of overhanging structures.
     
    756755                               THEN
    757756                                  pos_x = pos_x + pdx(i)
    758                                   CYCLE xloop                               
     757                                  CYCLE xloop
    759758                               ENDIF
    760759
     
    840839             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    841840
    842              DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
     841             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
    843842
    844843                particles(n)%id = 10000_idp**3 * grid_particles(kp,jp,ip)%id_counter + &
     
    857856!
    858857!-- Initialize aerosol background spectrum
    859     IF ( curvature_solution_effects  .AND.  .NOT. monodisperse_aerosols )  THEN
     858    IF ( curvature_solution_effects )  THEN
    860859       CALL lpm_init_aerosols(local_start)
    861860    ENDIF
     
    871870                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    872871!
    873 !--             Move only new particles. Moreover, limit random fluctuation 
    874 !--             in order to prevent that particles move more than one grid box, 
    875 !--             which would lead to problems concerning particle exchange 
    876 !--             between processors in case pdx/pdy are larger than dx/dy, 
    877 !--             respectively. 
     872!--             Move only new particles. Moreover, limit random fluctuation
     873!--             in order to prevent that particles move more than one grid box,
     874!--             which would lead to problems concerning particle exchange
     875!--             between processors in case pdx/pdy are larger than dx/dy,
     876!--             respectively.
    878877                DO  n = local_start(kp,jp,ip), number_of_particles
    879878                   IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
     
    883882                              MERGE( rand_contr, SIGN( dx, rand_contr ),       &
    884883                                     ABS( rand_contr ) < dx                    &
    885                                    ) 
     884                                   )
    886885                   ENDIF
    887886                   IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
     
    891890                              MERGE( rand_contr, SIGN( dy, rand_contr ),       &
    892891                                     ABS( rand_contr ) < dy                    &
    893                                    ) 
     892                                   )
    894893                   ENDIF
    895894                   IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
     
    899898                              MERGE( rand_contr, SIGN( dz, rand_contr ),       &
    900899                                     ABS( rand_contr ) < dz                    &
    901                                    ) 
     900                                   )
    902901                   ENDIF
    903902                ENDDO
     
    908907!
    909908!--             Furthermore, remove particles located in topography. Note, as
    910 !--             the particle speed is still zero at this point, wall 
     909!--             the particle speed is still zero at this point, wall
    911910!--             reflection boundary conditions will not work in this case.
    912911                particles =>                                                   &
     
    934933    ENDIF
    935934!
    936 !-- In case of random_start_position, delete particles identified by 
    937 !-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks, 
    938 !-- which is needed for a fast interpolation of the LES fields on the particle 
     935!-- In case of random_start_position, delete particles identified by
     936!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
     937!-- which is needed for a fast interpolation of the LES fields on the particle
    939938!-- position.
    940939    CALL lpm_pack_all_arrays
     
    967966
    968967    USE arrays_3d,                                                             &
    969         ONLY: hyp, pt, q 
     968        ONLY: hyp, pt, q
    970969
    971970    USE cloud_parameters,                                                      &
     
    978977
    979978    USE particle_attributes,                                                   &
    980         ONLY: init_aerosol_probabilistic, molecular_weight_of_solute,          &
    981               molecular_weight_of_water, n1, n2, n3, rho_s, rm1, rm2, rm3,     &
    982               s1, s2, s3, vanthoff
     979        ONLY: aero_type, aero_weight, log_sigma, molecular_weight_of_solute,   &
     980              molecular_weight_of_water, na, rho_s, rm, vanthoff
    983981
    984982    IMPLICIT NONE
    985 
    986     REAL(wp), DIMENSION(:), ALLOCATABLE ::  cdf     !< CDF of aerosol spectrum
    987     REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_temp  !< dry aerosol radius spectrum
    988983
    989984    REAL(wp)  :: afactor            !< curvature effects
    990985    REAL(wp)  :: bfactor            !< solute effects
    991     REAL(wp)  :: dr                 !< width of radius bin
     986    REAL(wp)  :: dlogr              !< logarithmic width of radius bin
    992987    REAL(wp)  :: e_a                !< vapor pressure
    993988    REAL(wp)  :: e_s                !< saturation vapor pressure
    994     REAL(wp)  :: n_init             !< sum of all aerosol concentrations
    995     REAL(wp)  :: pdf                !< PDF of aerosol spectrum
    996     REAL(wp)  :: rmin = 1.0e-8_wp   !< minimum aerosol radius
    997     REAL(wp)  :: rmax = 1.0e-6_wp   !< maximum aerosol radius
    998     REAL(wp)  :: rs_rand            !< random number
    999     REAL(wp)  :: r_mid              !< mean radius
     989    REAL(wp)  :: rmin = 0.005e-6_wp !< minimum aerosol radius
     990    REAL(wp)  :: rmax = 10.0e-6_wp  !< maximum aerosol radius
     991    REAL(wp)  :: r_mid              !< mean radius of bin
     992    REAL(wp)  :: r_l                !< left radius of bin
     993    REAL(wp)  :: r_r                !< right radius of bin
    1000994    REAL(wp)  :: sigma              !< surface tension
    1001995    REAL(wp)  :: t_int              !< temperature
    1002     REAL(wp)  :: weight_sum         !< sum of all weighting factors
    1003996
    1004997    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  local_start !<
    1005998
    1006999    INTEGER(iwp)  :: n              !<
    1007     INTEGER(iwp)  :: nn             !<
    1008     INTEGER(iwp)  :: no_bins = 999  !< number of bins
    10091000    INTEGER(iwp)  :: ip             !<
    10101001    INTEGER(iwp)  :: jp             !<
    10111002    INTEGER(iwp)  :: kp             !<
    10121003
    1013     LOGICAL ::  new_pdf = .FALSE.   !< check if aerosol PDF has to be recalculated
    1014 
    1015 !
    1016 !-- Compute aerosol background distribution
    1017     IF ( init_aerosol_probabilistic )  THEN
    1018        ALLOCATE( cdf(0:no_bins), r_temp(0:no_bins) )
    1019        DO n = 0, no_bins
    1020           r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) /            &
    1021                            REAL(no_bins, KIND=wp) * REAL(n, KIND=wp) )
    1022 
    1023           cdf(n) = 0.0_wp
    1024           n_init = n1 + n2 + n3
    1025           IF ( n1 > 0.0_wp )  THEN
    1026              cdf(n) = cdf(n) + n1 / n_init * ( 0.5_wp + 0.5_wp *        &
    1027                                   ERF( LOG( r_temp(n) / rm1 ) /         &
    1028                                        ( SQRT(2.0_wp) * LOG(s1) )       &
    1029                                      ) )
    1030           ENDIF
    1031           IF ( n2 > 0.0_wp )  THEN
    1032              cdf(n) = cdf(n) + n2 / n_init * ( 0.5_wp + 0.5_wp *        &
    1033                                   ERF( LOG( r_temp(n) / rm2 ) /         &
    1034                                        ( SQRT(2.0_wp) * LOG(s2) )       &
    1035                                      ) )
    1036           ENDIF
    1037           IF ( n3 > 0.0_wp )  THEN
    1038              cdf(n) = cdf(n) + n3 / n_init * ( 0.5_wp + 0.5_wp *        &
    1039                                   ERF( LOG( r_temp(n) / rm3 ) /         &
    1040                                        ( SQRT(2.0_wp) * LOG(s3) )       &
    1041                                      ) )
    1042           ENDIF
    1043 
    1044        ENDDO
     1004!
     1005!-- The following typical aerosol spectra are taken from Jaenicke (1993):
     1006!-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
     1007    IF ( TRIM(aero_type) .EQ. 'polar' )  THEN
     1008       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6
     1009       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6
     1010       log_sigma = (/ 0.245, 0.300, 0.291 /)
     1011    ELSEIF ( TRIM(aero_type) .EQ. 'background' )  THEN
     1012       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6
     1013       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6
     1014       log_sigma = (/ 0.645, 0.253, 0.425 /)
     1015    ELSEIF ( TRIM(aero_type) .EQ. 'maritime' )  THEN
     1016       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6
     1017       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6
     1018       log_sigma = (/ 0.657, 0.210, 0.396 /)
     1019    ELSEIF ( TRIM(aero_type) .EQ. 'continental' )  THEN
     1020       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6
     1021       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6
     1022       log_sigma = (/ 0.161, 0.217, 0.380 /)
     1023    ELSEIF ( TRIM(aero_type) .EQ. 'desert' )  THEN
     1024       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6
     1025       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6
     1026       log_sigma = (/ 0.247, 0.770, 0.438 /)
     1027    ELSEIF ( TRIM(aero_type) .EQ. 'rural' )  THEN
     1028       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6
     1029       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6
     1030       log_sigma = (/ 0.225, 0.557, 0.266 /)
     1031    ELSEIF ( TRIM(aero_type) .EQ. 'urban' )  THEN
     1032       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6
     1033       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6
     1034       log_sigma = (/ 0.245, 0.666, 0.337 /)
     1035    ELSEIF ( TRIM(aero_type) .EQ. 'user' )  THEN
     1036       CONTINUE
     1037    ELSE
     1038       WRITE( message_string, * ) 'unknown aerosol type ',   &
     1039                                'aero_type = "', TRIM( aero_type ), '"'
     1040       CALL message( 'lpm_init', 'PA0459', 1, 2, 0, 6, 0 )
    10451041    ENDIF
    10461042
     
    10521048             IF ( number_of_particles <= 0 )  CYCLE
    10531049             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     1050
     1051             dlogr   = ( LOG10(rmax) - LOG10(rmin) ) / ( number_of_particles - local_start(kp,jp,ip) + 1 )
    10541052!
    10551053!--          Initialize the aerosols with a predefined spectral distribution
    1056 !--          of the dry radius (logarithmically increasing bins) and a varying 
     1054!--          of the dry radius (logarithmically increasing bins) and a varying
    10571055!--          weighting factor
    1058              IF ( .NOT. init_aerosol_probabilistic )  THEN
    1059 
    1060                 new_pdf = .FALSE.
    1061                 IF ( .NOT. ALLOCATED( r_temp ) )  THEN
    1062                    new_pdf = .TRUE.
     1056             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
     1057
     1058                r_l   = 10.0**( LOG10( rmin ) + (n-1) * dlogr )
     1059                r_r   = 10.0**( LOG10( rmin ) + n * dlogr )
     1060                r_mid = SQRT( r_l * r_r )
     1061
     1062                particles(n)%aux1          = r_mid
     1063                particles(n)%weight_factor =                                           &
     1064                   ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) *                     &
     1065                     EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) +  &
     1066                     na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) *                     &
     1067                     EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) +  &
     1068                     na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) *                     &
     1069                     EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) )    &
     1070                   ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dz )
     1071
     1072!
     1073!--             Multiply weight_factor with the namelist parameter aero_weight
     1074!--             to increase or decrease the number of simulated aerosols
     1075                particles(n)%weight_factor = particles(n)%weight_factor * aero_weight
     1076
     1077                IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) &
     1078                     .GT. random_function( iran_part ) )  THEN
     1079                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0
    10631080                ELSE
    1064                    IF ( SIZE( r_temp ) .NE. &
    1065                         number_of_particles - local_start(kp,jp,ip) + 2 )  THEN
    1066                       new_pdf = .TRUE.
    1067                       DEALLOCATE( r_temp )
    1068                    ENDIF
     1081                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp)
    10691082                ENDIF
    1070 
    1071                 IF ( new_pdf )  THEN
    1072 
    1073                    no_bins = number_of_particles + 1 - local_start(kp,jp,ip)
    1074                    ALLOCATE( r_temp(0:no_bins) )
    1075 
    1076                    DO n = 0, no_bins
    1077                       r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / &
    1078                                        REAL(no_bins, KIND=wp) *                 &
    1079                                        REAL(n, KIND=wp) )
    1080                    ENDDO
    1081 
    1082                 ENDIF
    1083 
    1084 !
    1085 !--             Calculate radius and concentration of each aerosol
    1086                 DO n = local_start(kp,jp,ip), number_of_particles
    1087 
    1088                    nn = n - local_start(kp,jp,ip)
    1089 
    1090                    r_mid = SQRT( r_temp(nn) * r_temp(nn+1) )
    1091                    dr    = r_temp(nn+1) - r_temp(nn)
    1092 
    1093                    pdf    = 0.0_wp
    1094                    n_init = n1 + n2 + n3
    1095                    IF ( n1 > 0.0_wp )  THEN
    1096                       pdf = pdf + n1 / n_init * ( 1.0_wp / ( r_mid * LOG(s1) *      &
    1097                                                              SQRT( 2.0_wp * pi )    &
    1098                                                            ) *                      &
    1099                                                   EXP( -( LOG( r_mid / rm1 ) )**2 / &
    1100                                                        ( 2.0_wp * LOG(s1)**2 )      &
    1101                                                      )                              &
    1102                                                 )
    1103                    ENDIF
    1104                    IF ( n2 > 0.0_wp )  THEN
    1105                       pdf = pdf + n2 / n_init * ( 1.0_wp / ( r_mid * LOG(s2) *      &
    1106                                                              SQRT( 2.0_wp * pi )    &
    1107                                                            ) *                      &
    1108                                                   EXP( -( LOG( r_mid / rm2 ) )**2 / &
    1109                                                        ( 2.0_wp * LOG(s2)**2 )      &
    1110                                                      )                              &
    1111                                                 )
    1112                    ENDIF
    1113                    IF ( n3 > 0.0_wp )  THEN
    1114                       pdf = pdf + n3 / n_init * ( 1.0_wp / ( r_mid * LOG(s3) *      &
    1115                                                              SQRT( 2.0_wp * pi )    &
    1116                                                            ) *                      &
    1117                                                   EXP( -( LOG( r_mid / rm3 ) )**2 / &
    1118                                                        ( 2.0_wp * LOG(s3)**2 )      &
    1119                                                      )                              &
    1120                                                 )
    1121                    ENDIF
    1122 
    1123                    particles(n)%rvar2         = r_mid
    1124                    particles(n)%weight_factor = pdf * dr
    1125 
    1126                 END DO
    1127 !
    1128 !--             Adjust weighting factors to initialize the same number of aerosols
    1129 !--             in every grid box
    1130                 weight_sum = SUM(particles(local_start(kp,jp,ip):number_of_particles)%weight_factor)
    1131 
    1132                 particles(local_start(kp,jp,ip):number_of_particles)%weight_factor =     &
    1133                    particles(local_start(kp,jp,ip):number_of_particles)%weight_factor /  &
    1134                    weight_sum * initial_weighting_factor * ( no_bins + 1 )
    1135 
    1136              ENDIF
    1137 !
    1138 !--          Initialize the aerosols with a predefined weighting factor but
    1139 !--          a randomly choosen dry radius
    1140              IF ( init_aerosol_probabilistic )  THEN
    1141 
    1142                 DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
    1143 
    1144                    rs_rand = -1.0_wp
    1145                    DO WHILE ( rs_rand .LT. cdf(0)  .OR.  rs_rand .GE. cdf(no_bins)  )
    1146                       rs_rand = random_function( iran_part )
    1147                    ENDDO
    1148 !
    1149 !--                Determine aerosol dry radius by a random number generator
    1150                    DO nn = 0, no_bins-1
    1151                       IF ( cdf(nn) .LE. rs_rand  .AND.  cdf(nn+1) .GT. rs_rand )  THEN
    1152                          particles(n)%rvar2 = r_temp(nn) + ( r_temp(nn+1) - r_temp(nn) ) / &
    1153                                               ( cdf(nn+1) - cdf(nn) ) * ( rs_rand - cdf(nn) )
    1154                          EXIT
    1155                       ENDIF
    1156                    ENDDO
    1157 
    1158                 ENDDO
    1159 
    1160              ENDIF
    1161 
     1083!
     1084!--             Unnecessary particles will be deleted
     1085                IF ( particles(n)%weight_factor .LE. 0.0 )  particles(n)%particle_mask = .FALSE.
     1086
     1087             ENDDO
    11621088!
    11631089!--          Set particle radius to equilibrium radius based on the environmental
    11641090!--          supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids
    1165 !--          the sometimes lengthy growth toward their equilibrium radius within 
     1091!--          the sometimes lengthy growth toward their equilibrium radius within
    11661092!--          the simulation.
    11671093             t_int  = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
     
    11761102                       rho_s / ( molecular_weight_of_solute * rho_l )
    11771103!
    1178 !--          The formula is only valid for subsaturated environments. For 
     1104!--          The formula is only valid for subsaturated environments. For
    11791105!--          supersaturations higher than -5 %, the supersaturation is set to -5%.
    11801106             IF ( e_a / e_s >= 0.95_wp )  e_a = 0.95_wp * e_s
    11811107
    1182              DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
    1183 !
    1184 !--             For details on this equation, see Eq. (14) of Khvorostyanov and 
    1185 !--             Curry (2007, JGR) 
     1108             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
     1109!
     1110!--             For details on this equation, see Eq. (14) of Khvorostyanov and
     1111!--             Curry (2007, JGR)
    11861112                particles(n)%radius = bfactor**0.3333333_wp *                  &
    1187                    particles(n)%rvar2 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
     1113                   particles(n)%aux1 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
    11881114                   ( 1.0_wp + ( afactor / ( 3.0_wp * bfactor**0.3333333_wp *   &
    1189                      particles(n)%rvar2 ) ) /                                  &
     1115                     particles(n)%aux1 ) ) /                                  &
    11901116                     ( 1.0_wp - e_a / e_s )**0.6666666_wp                      &
    11911117                   )
     
    11961122       ENDDO
    11971123    ENDDO
    1198 !
    1199 !-- Deallocate used arrays
    1200     IF ( ALLOCATED(r_temp) )  DEALLOCATE( r_temp )
    1201     IF ( ALLOCATED(cdf) )     DEALLOCATE( cdf )
    12021124
    12031125 END SUBROUTINE lpm_init_aerosols
  • palm/trunk/SOURCE/lpm_read_restart_file.f90

    r2305 r2312  
    2020! Current revisions:
    2121! ------------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Extended particle data type.
     28!
     29! 2305 2017-07-06 11:18:47Z hoffmann
    2730! Improved calculation of particle IDs.
    28 ! 
     31!
    2932! 2265 2017-06-08 16:58:28Z schwenkel
    3033! Unused variables removed.
    31 ! 
     34!
    3235! 2101 2017-01-05 16:42:31Z suehring
    3336!
    3437! 2000 2016-08-20 18:09:15Z knoop
    3538! Forced header and separation lines into 80 columns
    36 ! 
     39!
    3740! 1822 2016-04-07 07:49:42Z hoffmann
    3841! Tails removed. Unused variables removed.
    3942!
    4043! 1682 2015-10-07 23:56:08Z knoop
    41 ! Code annotations made doxygen readable 
    42 ! 
     44! Code annotations made doxygen readable
     45!
    4346! 1359 2014-04-11 17:15:14Z hoffmann
    44 ! New particle structure integrated. 
     47! New particle structure integrated.
    4548!
    4649! 1320 2014-03-20 08:40:49Z raasch
     
    6164!------------------------------------------------------------------------------!
    6265 SUBROUTINE lpm_read_restart_file
    63  
     66
    6467
    6568    USE control_parameters,                                                    &
     
    119122
    120123!
    121 !-- If less particles are stored on the restart file than prescribed by 
     124!-- If less particles are stored on the restart file than prescribed by
    122125!-- min_nr_particle, the remainder is initialized by zero_particle to avoid
    123126!-- errors.
     
    125128                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    126129                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    127                                    0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
     130                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    128131                                   0, 0, 0_idp, .FALSE., -1 )
    129132!
     
    173176    CLOSE ( 90 )
    174177!
    175 !-- Must be called to sort particles into blocks, which is needed for a fast 
     178!-- Must be called to sort particles into blocks, which is needed for a fast
    176179!-- interpolation of the LES fields on the particle position.
    177180    CALL lpm_pack_all_arrays
  • palm/trunk/SOURCE/microphysics_mod.f90

    r2292 r2312  
    2020! Current revisions:
    2121! ------------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
    27 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
    28 ! includes two more prognostic equations for cloud drop concentration (nc) 
    29 ! and cloud water content (qc).
    30 ! - The process of activation is parameterized with a simple Twomey
    31 !   activion scheme or with considering solution and curvature
     27! s1 changed to log_sigma
     28!
     29! 2292 2017-06-20 09:51:42Z schwenkel
     30! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     31! includes two more prognostic equations for cloud drop concentration (nc)
     32! and cloud water content (qc).
     33! - The process of activation is parameterized with a simple Twomey
     34!   activion scheme or with considering solution and curvature
    3235!   effects (Khvorostyanov and Curry ,2006).
    3336! - The saturation adjustment scheme is replaced by the parameterization
    3437!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
    3538! - All other microphysical processes of Seifert and Beheng are used.
    36 !   Additionally, in those processes the reduction of cloud number concentration 
    37 !   is considered. 
    38 ! 
     39!   Additionally, in those processes the reduction of cloud number concentration
     40!   is considered.
     41!
    3942! 2233 2017-05-30 18:08:54Z suehring
    4043!
    4144! 2232 2017-05-30 17:47:52Z suehring
    4245! Adjustments to new topography and surface concept
    43 ! 
     46!
    4447! 2155 2017-02-21 09:57:40Z hoffmann
    4548! Bugfix in the calculation of microphysical quantities on ghost points.
     
    207210           limiter_sedimentation, na_init, nc_const, sigma_gc,                 &
    208211           ventilation_effect
    209            
     212
    210213
    211214    INTERFACE microphysics_control
     
    463466                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
    464467                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin *               &
    465                                        MERGE( 1.0_wp, 0.0_wp,                  &           
     468                                       MERGE( 1.0_wp, 0.0_wp,                  &
    466469                                              BTEST( wall_flags_0(k,j,i), 0 ) )
    467470                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
    468471                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax *               &
    469                                        MERGE( 1.0_wp, 0.0_wp,                  &           
     472                                       MERGE( 1.0_wp, 0.0_wp,                  &
    470473                                              BTEST( wall_flags_0(k,j,i), 0 ) )
    471474                   ENDIF
     
    478481                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
    479482                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin *            &
    480                                           MERGE( 1.0_wp, 0.0_wp,               &           
     483                                          MERGE( 1.0_wp, 0.0_wp,               &
    481484                                              BTEST( wall_flags_0(k,j,i), 0 ) )
    482485                      ENDIF
     
    521524       USE particle_attributes,                                                &
    522525          ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
    523               s1, s2, s3, vanthoff
     526              log_sigma, vanthoff
    524527
    525528       IMPLICIT NONE
     
    578581!--             Prescribe parameters for activation
    579582!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
    580                 k_act  = 0.7_wp 
     583                k_act  = 0.7_wp
    581584
    582585                IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
    583586!
    584 !--                Compute the number of activated Aerosols 
     587!--                Compute the number of activated Aerosols
    585588!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
    586589                   n_act     = na_init * sat**k_act
     
    592595!
    593596!--                Compute activation rate after Khairoutdinov and Kogan
    594 !--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
    595                    sat_max = 1.0_wp / 100.0_wp                 
    596                    activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       & 
    597                       ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    & 
     597!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
     598                   sat_max = 1.0_wp / 100.0_wp
     599                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
     600                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
    598601                       dt_micro
    599602                ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
    600603!
    601 !--                Curvature effect (afactor) with surface tension 
     604!--                Curvature effect (afactor) with surface tension
    602605!--                parameterization by Straka (2009)
    603606                   sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
     
    609612
    610613!
    611 !--                Prescribe power index that describes the soluble fraction 
     614!--                Prescribe power index that describes the soluble fraction
    612615!--                of an aerosol particle (beta) and mean geometric radius of
    613616!--                dry aerosol spectrum
     
    616619                   rd0      = 0.05E-6_wp
    617620!
    618 !--                Calculate mean geometric supersaturation (s_0) with 
     621!--                Calculate mean geometric supersaturation (s_0) with
    619622!--                parameterization by Khvorostyanov and Curry (2006)
    620623                   s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                  &
     
    622625
    623626!
    624 !--                Calculate number of activated CCN as a function of 
     627!--                Calculate number of activated CCN as a function of
    625628!--                supersaturation and taking Koehler theory into account
    626 !--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)       
    627                    n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              & 
    628                       LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )     
     629!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
     630                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
     631                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * log_sigma(1) ) ) )
    629632                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
    630633                ENDIF
     
    644647! Description:
    645648! ------------
    646 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 
     649!> Calculate condensation rate for cloud water content (after Khairoutdinov and
    647650!> Kogan, 2000).
    648651!------------------------------------------------------------------------------!
     
    724727!--             Mean weight of cloud drops
    725728                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
    726                 xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)               
     729                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
    727730!
    728731!--             Weight averaged diameter of cloud drops:
     
    730733!
    731734!--             Integral diameter of cloud drops
    732                 nc_0 = nc(k,j,i) * dc               
     735                nc_0 = nc(k,j,i) * dc
    733736!
    734737!--             Condensation needs only to be calculated in supersaturated regions
     
    741744                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
    742745                   cond      = MIN( cond, cond_max / dt_micro )
    743                
     746
    744747                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro
    745748                ELSEIF ( sat < 0.0_wp ) THEN
     
    891894                                                              * flag
    892895                   IF ( microphysics_morrison ) THEN
    893                       nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         & 
     896                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
    894897                                    autocon / x0 * hyrho(k) * dt_micro * flag )
    895898                   ENDIF
     
    10331036!
    10341037!--                Mean weight of cloud drops
    1035                    xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)     
     1038                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
    10361039!
    10371040!--                Parameterized turbulence effects on autoconversion (Seifert,
     
    13211324
    13221325       USE control_parameters,                                                 &
    1323            ONLY:  call_microphysics_at_all_substeps,                           & 
     1326           ONLY:  call_microphysics_at_all_substeps,                           &
    13241327                  intermediate_timestep_count, microphysics_morrison
    13251328
     
    13661369                IF ( microphysics_morrison ) THEN
    13671370                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
    1368                       sed_nc(k) = sed_qc_const *                               & 
     1371                      sed_nc(k) = sed_qc_const *                               &
    13691372                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
    13701373                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
     
    14491452       USE surface_mod,                                                        &
    14501453           ONLY :  bc_h
    1451        
     1454
    14521455       IMPLICIT NONE
    14531456
     
    14551458       INTEGER(iwp) ::  j             !< running index y direction
    14561459       INTEGER(iwp) ::  k             !< running index z direction
    1457        INTEGER(iwp) ::  k_run         !< 
     1460       INTEGER(iwp) ::  k_run         !<
    14581461       INTEGER(iwp) ::  l             !< running index of surface type
    14591462       INTEGER(iwp) ::  m             !< running index surface elements
     
    15281531             ENDDO
    15291532!
    1530 !--          Adjust boundary values using surface data type. 
     1533!--          Adjust boundary values using surface data type.
    15311534!--          Upward-facing
    1532              surf_s = bc_h(0)%start_index(j,i)   
     1535             surf_s = bc_h(0)%start_index(j,i)
    15331536             surf_e = bc_h(0)%end_index(j,i)
    15341537             DO  m = surf_s, surf_e
     
    15391542!
    15401543!--          Downward-facing
    1541              surf_s = bc_h(1)%start_index(j,i)   
     1544             surf_s = bc_h(1)%start_index(j,i)
    15421545             surf_e = bc_h(1)%end_index(j,i)
    15431546             DO  m = surf_s, surf_e
     
    16071610                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    16081611!
    1609 !--             Sum up all rain drop number densities which contribute to the flux 
     1612!--             Sum up all rain drop number densities which contribute to the flux
    16101613!--             through k-1/2
    16111614                flux  = 0.0_wp
     
    17301733       THEN
    17311734!
    1732 !--       Run over all upward-facing surface elements, i.e. non-natural, 
     1735!--       Run over all upward-facing surface elements, i.e. non-natural,
    17331736!--       natural and urban
    17341737          DO  m = 1, bc_h(0)%ns
    1735              i = bc_h(0)%i(m)           
     1738             i = bc_h(0)%i(m)
    17361739             j = bc_h(0)%j(m)
    17371740             k = bc_h(0)%k(m)
     
    19611964       USE particle_attributes,                                                &
    19621965          ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
    1963               s1, s2, s3, vanthoff
     1966              log_sigma, vanthoff
    19641967
    19651968       IMPLICIT NONE
     
    20132016!--       Prescribe parameters for activation
    20142017!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
    2015           k_act  = 0.7_wp 
     2018          k_act  = 0.7_wp
    20162019
    20172020          IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
    20182021!
    2019 !--          Compute the number of activated Aerosols 
     2022!--          Compute the number of activated Aerosols
    20202023!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
    20212024             n_act     = na_init * sat**k_act
     
    20272030!
    20282031!--          Compute activation rate after Khairoutdinov and Kogan
    2029 !--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
    2030              sat_max = 0.8_wp / 100.0_wp                 
    2031              activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              & 
    2032                  ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          & 
     2032!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
     2033             sat_max = 0.8_wp / 100.0_wp
     2034             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              &
     2035                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          &
    20332036                  dt_micro
    20342037
     
    20362039          ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk )  THEN
    20372040!
    2038 !--          Curvature effect (afactor) with surface tension 
     2041!--          Curvature effect (afactor) with surface tension
    20392042!--          parameterization by Straka (2009)
    20402043             sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
     
    20462049
    20472050!
    2048 !--          Prescribe power index that describes the soluble fraction 
     2051!--          Prescribe power index that describes the soluble fraction
    20492052!--          of an aerosol particle (beta) and mean geometric radius of
    20502053!--          dry aerosol spectrum
     
    20532056             rd0      = 0.05E-6_wp
    20542057!
    2055 !--          Calculate mean geometric supersaturation (s_0) with 
     2058!--          Calculate mean geometric supersaturation (s_0) with
    20562059!--          parameterization by Khvorostyanov and Curry (2006)
    20572060             s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                        &
     
    20592062
    20602063!
    2061 !--          Calculate number of activated CCN as a function of 
     2064!--          Calculate number of activated CCN as a function of
    20622065!--          supersaturation and taking Koehler theory into account
    20632066!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
    2064              n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                     & 
    2065                 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )     
     2067             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                     &
     2068                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * log_sigma(1) ) ) )
    20662069             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
    2067                                      
    2068              nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)                   
     2070
     2071             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
    20692072          ENDIF
    20702073
     
    20762079! Description:
    20772080! ------------
    2078 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 
     2081!> Calculate condensation rate for cloud water content (after Khairoutdinov and
    20792082!> Kogan, 2000).
    20802083!------------------------------------------------------------------------------!
     
    21542157!--       Mean weight of cloud drops
    21552158          IF ( nc_1d(k) <= 0.0_wp) CYCLE
    2156           xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)               
     2159          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
    21572160!
    21582161!--       Weight averaged diameter of cloud drops:
     
    21602163!
    21612164!--       Integral diameter of cloud drops
    2162           nc_0 = nc_1d(k) * dc               
     2165          nc_0 = nc_1d(k) * dc
    21632166!
    21642167!--       Condensation needs only to be calculated in supersaturated regions
     
    21712174             cond_max  = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
    21722175             cond      = MIN( cond, cond_max / dt_micro )
    2173                
     2176
    21742177             qc_1d(k) = qc_1d(k) + cond * dt_micro
    21752178          ELSEIF ( sat < 0.0_wp ) THEN
     
    23022305             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
    23032306             IF ( microphysics_morrison ) THEN
    2304                 nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                & 
     2307                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                &
    23052308                                    autocon / x0 * hyrho(k) * dt_micro * flag )
    23062309             ENDIF
     
    24252428!
    24262429!--          Mean weight of cloud drops
    2427              xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)     
     2430             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
    24282431!
    24292432!--          Parameterized turbulence effects on autoconversion (Seifert,
     
    27162719          IF ( microphysics_morrison ) THEN
    27172720             IF ( qc_1d(k) > eps_sb  .AND.  nc_1d(k) > eps_mr )  THEN
    2718                 sed_nc(k) = sed_qc_const *                                     & 
     2721                sed_nc(k) = sed_qc_const *                                     &
    27192722                            ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
    27202723                            ( nc_1d(k) )**( 1.0_wp / 3.0_wp )
     
    27902793       USE surface_mod,                                                        &
    27912794           ONLY :  bc_h
    2792        
     2795
    27932796       IMPLICIT NONE
    27942797
     
    28592862       ENDDO
    28602863!
    2861 !--    Adjust boundary values using surface data type. 
     2864!--    Adjust boundary values using surface data type.
    28622865!--    Upward facing non-natural
    2863        surf_s = bc_h(0)%start_index(j,i)   
     2866       surf_s = bc_h(0)%start_index(j,i)
    28642867       surf_e = bc_h(0)%end_index(j,i)
    28652868       DO  m = surf_s, surf_e
     
    28702873!
    28712874!--    Downward facing non-natural
    2872        surf_s = bc_h(1)%start_index(j,i)   
     2875       surf_s = bc_h(1)%start_index(j,i)
    28732876       surf_e = bc_h(1)%end_index(j,i)
    28742877       DO  m = surf_s, surf_e
     
    29352938          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    29362939!
    2937 !--       Sum up all rain drop number densities which contribute to the flux 
     2940!--       Sum up all rain drop number densities which contribute to the flux
    29382941!--       through k-1/2
    29392942          flux  = 0.0_wp
     
    30443047       THEN
    30453048
    3046           surf_s = bc_h(0)%start_index(j,i)   
     3049          surf_s = bc_h(0)%start_index(j,i)
    30473050          surf_e = bc_h(0)%end_index(j,i)
    30483051          DO  m = surf_s, surf_e
  • palm/trunk/SOURCE/mod_particle_attributes.f90

    r2305 r2312  
    2020! Current revisions:
    2121! ------------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Aerosol initialization improved.
     28!
     29! 2305 2017-07-06 11:18:47Z hoffmann
    2730! Improved calculation of particle IDs.
    28 ! 
     31!
    2932! 2278 2017-06-12 13:08:18Z schwenkel
    3033! Added comments
    31 ! 
     34!
    3235! 2265 2017-06-08 16:58:28Z schwenkel
    3336! Unused variables removed.
    34 ! 
     37!
    3538! 2263 2017-06-08 14:59:01Z schwenkel
    3639! Implemented splitting and merging algorithm
    37 ! 
     40!
    3841! 2183 2017-03-17 14:29:15Z schwenkel
    3942!
    4043! 2182 2017-03-17 14:27:40Z schwenkel
    4144! Added parameters for simplified particle initialization.
    42 ! 
     45!
    4346! 2122 2017-01-18 12:22:54Z hoffmann
    4447! Calculation of particle ID
     
    4851! 2000 2016-08-20 18:09:15Z knoop
    4952! Forced header and separation lines into 80 columns
    50 ! 
     53!
    5154! 1936 2016-06-13 13:37:44Z suehring
    5255! +deallocate_memory, step_dealloc
     
    7174!
    7275! 1727 2015-11-20 07:22:02Z knoop
    73 ! Bugfix: Cause of syntax warning gfortran preprocessor removed 
    74 ! 
     76! Bugfix: Cause of syntax warning gfortran preprocessor removed
     77!
    7578! 1682 2015-10-07 23:56:08Z knoop
    76 ! Code annotations made doxygen readable 
     79! Code annotations made doxygen readable
    7780!
    7881! 1575 2015-03-27 09:56:27Z raasch
     
    8891!------------------------------------------------------------------------------!
    8992MODULE particle_attributes
    90  
     93
    9194
    9295    USE kinds
    9396
    94     CHARACTER(LEN=15) ::  bc_par_lr = 'cyclic'                    !< left/right boundary condition
    95     CHARACTER(LEN=15) ::  bc_par_ns = 'cyclic'                    !< north/south boundary condition
    96     CHARACTER(LEN=15) ::  bc_par_b  = 'reflect'                   !< bottom boundary condition
    97     CHARACTER(LEN=15) ::  bc_par_t  = 'absorb'                    !< top boundary condition
    98     CHARACTER(LEN=15) ::  collision_algorithm = 'all_or_nothing'  !< collision algorithm
    99     CHARACTER(LEN=15) ::  collision_kernel = 'none'               !< collision kernel
    100     CHARACTER(LEN=5) ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
    101     CHARACTER(LEN=5) ::  splitting_mode = 'const'                 !< splitting mode
    102 
    103 
    104     INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step
     97    CHARACTER(LEN=15) ::  aero_type = 'maritime'                   !< aerosol type
     98    CHARACTER(LEN=15) ::  bc_par_lr = 'cyclic'                     !< left/right boundary condition
     99    CHARACTER(LEN=15) ::  bc_par_ns = 'cyclic'                     !< north/south boundary condition
     100    CHARACTER(LEN=15) ::  bc_par_b  = 'reflect'                    !< bottom boundary condition
     101    CHARACTER(LEN=15) ::  bc_par_t  = 'absorb'                     !< top boundary condition
     102    CHARACTER(LEN=15) ::  collision_kernel = 'none'                !< collision kernel
     103    CHARACTER(LEN=5)  ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
     104    CHARACTER(LEN=5)  ::  splitting_mode = 'const'                 !< splitting mode
     105
     106    INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step
    105107    INTEGER(iwp) ::  dissipation_classes = 10                     !< namelist parameter (see documentation)
    106108    INTEGER(iwp) ::  ibc_par_b                                    !< particle bottom boundary condition dummy
     
    113115    INTEGER(iwp) ::  max_number_particles_per_gridbox = 100       !< namelist parameter (see documentation)
    114116    INTEGER(iwp) ::  merge_drp = 0                                !< number of merged droplets
    115     INTEGER(iwp) ::  min_nr_particle = 50                         !< namelist parameter (see documentation)         
     117    INTEGER(iwp) ::  min_nr_particle = 50                         !< namelist parameter (see documentation)
    116118    INTEGER(iwp) ::  new_particles = 0                            !< number of new particles
    117     INTEGER(iwp) ::  n_max = 100                                  !< number of radii bin for splitting functions     
    118     INTEGER(iwp) ::  number_of_particles = 0                      !< number of particles for each grid box (3d array is saved on prt_count)           
    119     INTEGER(iwp) ::  number_of_particle_groups = 1                !< namelist parameter (see documentation)           
    120     INTEGER(iwp) ::  number_of_sublayers = 20                     !< number of sublayers for particle velocities betwenn surface and first grid level         
    121     INTEGER(iwp) ::  number_particles_per_gridbox = -1            !< namelist parameter (see documentation)         
    122     INTEGER(iwp) ::  offset_ocean_nzt = 0                         !< in case of oceans runs, the vertical index calculations need an offset         
     119    INTEGER(iwp) ::  n_max = 100                                  !< number of radii bin for splitting functions
     120    INTEGER(iwp) ::  number_of_particles = 0                      !< number of particles for each grid box (3d array is saved on prt_count)
     121    INTEGER(iwp) ::  number_of_particle_groups = 1                !< namelist parameter (see documentation)
     122    INTEGER(iwp) ::  number_of_sublayers = 20                     !< number of sublayers for particle velocities betwenn surface and first grid level
     123    INTEGER(iwp) ::  number_particles_per_gridbox = -1            !< namelist parameter (see documentation)
     124    INTEGER(iwp) ::  offset_ocean_nzt = 0                         !< in case of oceans runs, the vertical index calculations need an offset
    123125    INTEGER(iwp) ::  offset_ocean_nzt_m1 = 0                      !< in case of oceans runs, the vertical index calculations need an offset
    124126    INTEGER(iwp) ::  particles_per_point = 1                      !< namelist parameter (see documentation)
     
    130132    INTEGER(iwp) ::  sum_merge_drp = 0                            !< sum of merged super droplets
    131133    INTEGER(iwp) ::  sum_new_particles = 0                        !< sum of created particles (in splitting algorithm)
    132     INTEGER(iwp) ::  total_number_of_particles                    !< total number of particles in the whole model domain           
     134    INTEGER(iwp) ::  total_number_of_particles                    !< total number of particles in the whole model domain
    133135    INTEGER(iwp) ::  trlp_count_sum                               !< parameter for particle exchange of PEs
    134136    INTEGER(iwp) ::  trlp_count_recv_sum                          !< parameter for particle exchange of PEs
     
    144146    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  prt_count  !< 3d array of number of particles of every grid box
    145147
    146     LOGICAL ::  all_or_nothing = .FALSE.                  !< flag for collision algorithm
    147     LOGICAL ::  average_impact = .FALSE.                  !< flag for collision algortihm
    148     LOGICAL ::  curvature_solution_effects = .FALSE.      !< namelist parameter (see documentation)                     
    149     LOGICAL ::  deallocate_memory = .TRUE.                !< namelist parameter (see documentation)                 
     148    LOGICAL ::  curvature_solution_effects = .FALSE.      !< namelist parameter (see documentation)
     149    LOGICAL ::  deallocate_memory = .TRUE.                !< namelist parameter (see documentation)
    150150    LOGICAL ::  hall_kernel = .FALSE.                     !< flag for collision kernel
    151     LOGICAL ::  init_aerosol_probabilistic = .FALSE.      !< namelist parameter (see documentation)
    152151    LOGICAL ::  merging = .FALSE.                         !< namelist parameter (see documentation)
    153     LOGICAL ::  monodisperse_aerosols = .FALSE.           !< namelist parameter (see documentation)
    154152    LOGICAL ::  particle_advection = .FALSE.              !< parameter to steer the advection of particles
    155     LOGICAL ::  random_start_position = .FALSE.           !< namelist parameter (see documentation)                   
    156     LOGICAL ::  read_particles_from_restartfile = .TRUE.  !< namelist parameter (see documentation)                   
     153    LOGICAL ::  random_start_position = .FALSE.           !< namelist parameter (see documentation)
     154    LOGICAL ::  read_particles_from_restartfile = .TRUE.  !< namelist parameter (see documentation)
    157155    LOGICAL ::  seed_follows_topography = .FALSE.         !< namelist parameter (see documentation)
    158156    LOGICAL ::  splitting = .FALSE.                       !< namelist parameter (see documentation)
    159157    LOGICAL ::  use_kernel_tables = .FALSE.               !< parameter, which turns on the use of precalculated collision kernels
    160     LOGICAL ::  use_sgs_for_particles = .FALSE.           !< namelist parameter (see documentation) 
     158    LOGICAL ::  use_sgs_for_particles = .FALSE.           !< namelist parameter (see documentation)
    161159    LOGICAL ::  wang_kernel = .FALSE.                     !< flag for collision kernel
    162160    LOGICAL ::  write_particle_statistics = .FALSE.       !< namelist parameter (see documentation)
    163                
     161
    164162    LOGICAL, DIMENSION(max_number_of_particle_groups) ::                       &
    165163                vertical_particle_advection = .TRUE.              !< Switch on/off vertical particle transport
    166164
     165    REAL(wp) ::  aero_weight = 1.0_wp                      !< namelist parameter (see documentation)
    167166    REAL(wp) ::  alloc_factor = 20.0_wp                    !< namelist parameter (see documentation)
    168167    REAL(wp) ::  c_0 = 3.0_wp                              !< parameter for lagrangian timescale
    169168    REAL(wp) ::  dt_min_part = 0.0002_wp                   !< minimum particle time step when SGS velocities are used (s)
    170169    REAL(wp) ::  dt_prel = 9999999.9_wp                    !< namelist parameter (see documentation)
    171     REAL(wp) ::  dt_write_particle_data = 9999999.9_wp     !< namelist parameter (see documentation)           
    172     REAL(wp) ::  end_time_prel = 9999999.9_wp              !< namelist parameter (see documentation)         
     170    REAL(wp) ::  dt_write_particle_data = 9999999.9_wp     !< namelist parameter (see documentation)
     171    REAL(wp) ::  end_time_prel = 9999999.9_wp              !< namelist parameter (see documentation)
    173172    REAL(wp) ::  initial_weighting_factor = 1.0_wp         !< namelist parameter (see documentation)
     173    REAL(wp) ::  log_sigma(3) = 1.0_wp                     !< namelist parameter (see documentation)
    174174    REAL(wp) ::  molecular_weight_of_solute = 0.05844_wp   !< mol. m. NaCl (kg mol-1)
    175175    REAL(wp) ::  molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1)
    176     REAL(wp) ::  n1 = 100.0_wp                             !< namelist parameter (see documentation)
    177     REAL(wp) ::  n2 = 0.0_wp                               !< namelist parameter (see documentation)
    178     REAL(wp) ::  n3 = 0.0_wp                               !< namelist parameter (see documentation)
     176    REAL(wp) ::  na(3) = 0.0_wp                            !< namelist parameter (see documentation)
    179177    REAL(wp) ::  number_concentration = -1.0_wp            !< namelist parameter (see documentation)
    180178    REAL(wp) ::  particle_advection_start = 0.0_wp         !< namelist parameter (see documentation)
    181179    REAL(wp) ::  radius_merge = 1.0E-7_wp                  !< namelist parameter (see documentation)
    182180    REAL(wp) ::  radius_split = 40.0E-6_wp                 !< namelist parameter (see documentation)
    183     REAL(wp) ::  rho_s = 2165.0_wp                         !< density of NaCl (kg m-3)
    184     REAL(wp) ::  rm1 = 0.05E-6_wp                          !< namelist parameter (see documentation)
    185     REAL(wp) ::  rm2 = 0.05E-6_wp                          !< namelist parameter (see documentation)
    186     REAL(wp) ::  rm3 = 0.05E-6_wp                          !< namelist parameter (see documentation)
    187     REAL(wp) ::  s1 = 2.0_wp                               !< namelist parameter (see documentation)
    188     REAL(wp) ::  s2 = 2.0_wp                               !< namelist parameter (see documentation)
    189     REAL(wp) ::  s3 = 2.0_wp                               !< namelist parameter (see documentation)
    190     REAL(wp) ::  sgs_wf_part                               !< parameter for sgs
     181    REAL(wp) ::  rho_s = 2165.0_wp                         !< density of NaCl (kg m-3)
     182    REAL(wp) ::  rm(3) = 1.0E-6_wp                         !< namelist parameter (see documentation)
     183    REAL(wp) ::  sgs_wf_part                               !< parameter for sgs
    191184    REAL(wp) ::  time_prel = 0.0_wp                        !< time for particle release
    192185    REAL(wp) ::  time_write_particle_data = 0.0_wp         !< write particle data at current time on file
     
    210203    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0)
    211204
    212    
     205
    213206    TYPE particle_type
     207        REAL(wp)     ::  aux1          !< auxiliary multi-purpose feature
     208        REAL(wp)     ::  aux2          !< auxiliary multi-purpose feature
    214209        REAL(wp)     ::  radius        !< radius of particle
    215210        REAL(wp)     ::  age           !< age of particle
    216         REAL(wp)     ::  age_m         !< 
     211        REAL(wp)     ::  age_m         !<
    217212        REAL(wp)     ::  dt_sum        !<
    218         REAL(wp)     ::  user          !< varible for user
    219         REAL(wp)     ::  e_m           !< interpolated sgs tke     
     213        REAL(wp)     ::  e_m           !< interpolated sgs tke
    220214        REAL(wp)     ::  origin_x      !< origin x-position of particle (changed cyclic bc)
    221215        REAL(wp)     ::  origin_y      !< origin y-position of particle (changed cyclic bc)
    222216        REAL(wp)     ::  origin_z      !< origin z-position of particle (changed cyclic bc)
    223         REAL(wp)     ::  rvar1         !< 
     217        REAL(wp)     ::  rvar1         !<
    224218        REAL(wp)     ::  rvar2         !<
    225         REAL(wp)     ::  rvar3         !< 
     219        REAL(wp)     ::  rvar3         !<
    226220        REAL(wp)     ::  speed_x       !< speed of particle in x
    227221        REAL(wp)     ::  speed_y       !< speed of particle in y
    228222        REAL(wp)     ::  speed_z       !< speed of particle in z
    229         REAL(wp)     ::  weight_factor !< weighting factor 
    230         REAL(wp)     ::  x             !< x-position 
    231         REAL(wp)     ::  y             !< y-position 
    232         REAL(wp)     ::  z             !< z-position 
     223        REAL(wp)     ::  weight_factor !< weighting factor
     224        REAL(wp)     ::  x             !< x-position
     225        REAL(wp)     ::  y             !< y-position
     226        REAL(wp)     ::  z             !< z-position
    233227        INTEGER(iwp) ::  class         !< radius class needed for collision
    234228        INTEGER(iwp) ::  group         !< number of particle group
     
    274268
    275269END MODULE particle_attributes
    276 
  • palm/trunk/SOURCE/package_parin.f90

    r2263 r2312  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Aerosol initialization improved.
     28!
     29! 2263 2017-06-08 14:59:01Z schwenkel
    2730! Implemented splitting and merging algorithm
    28 ! 
     31!
    2932! 2183 2017-03-17 14:29:15Z schwenkel
    3033!
    3134! 2182 2017-03-17 14:27:40Z schwenkel
    3235! Added parameters for simplified particle initialization.
    33 ! 
     36!
    3437! 2000 2016-08-20 18:09:15Z knoop
    3538! Forced header and separation lines into 80 columns
    36 ! 
     39!
    3740! 1936 2016-06-13 13:37:44Z suehring
    3841! +deallocate_memory, step_dealloc
    39 ! 
     42!
    4043! 1871 2016-04-15 11:46:09Z hoffmann
    4144! Initialization of aerosols added.
    42 ! 
     45!
    4346! 1833 2016-04-07 14:23:03Z raasch
    4447! reading of spectra_par moved to spectra_mod
    45 ! 
     48!
    4649! 1831 2016-04-07 13:15:51Z hoffmann
    4750! curvature_solution_effects added
     
    5053! Reading of &radiation_par moved to radiation_model_mod.
    5154! Reading of &canopy_par moved to plant_canopy_model_mod.
    52 ! 
     55!
    5356! 822 2016-04-07 07:49:42Z hoffmann
    5457! +collision_algorithm
     
    5760! 1817 2016-04-06 15:44:20Z maronga
    5861! Reading of &lsm_par moved to land_surface_model_mod.
    59 ! 
     62!
    6063! 1788 2016-03-10 11:01:04Z maronga
    6164! Parameter dewfall removed.
     
    6669! 1757 2016-02-22 15:49:32Z maronga
    6770! Added parameter unscheduled_radiation_calls
    68 ! 
     71!
    6972! 1691 2015-10-26 16:17:44Z maronga
    7073! Added skip_time_do_lsm, skip_time_do_radiation, and emissivity
    71 ! 
     74!
    7275! 1682 2015-10-07 23:56:08Z knoop
    73 ! Code annotations made doxygen readable 
    74 ! 
     76! Code annotations made doxygen readable
     77!
    7578! 1585 2015-04-30 07:05:52Z maronga
    7679! Added several radiation_par parameters
     
    8184! 1553 2015-03-03 17:33:54Z maronga
    8285! Resorting of lsm_par
    83 ! 
     86!
    8487! 1551 2015-03-03 14:18:16Z maronga
    8588! Several changes in the liste for land surface model and radiation model
    86 ! 
     89!
    8790! 1496 2014-12-02 17:25:50Z maronga
    8891! Added support for the land surface model and radiation scheme
    89 ! 
     92!
    9093! 1484 2014-10-21 10:53:05Z kanani
    9194! Changes due to new module structure of the plant canopy model:
    9295!   module plant_canopy_model_mod added,
    93 !   new package/namelist canopy_par added, i.e. the canopy model is no longer 
     96!   new package/namelist canopy_par added, i.e. the canopy model is no longer
    9497!   steered over the inipar-namelist,
    9598!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
    9699!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff.
    97100! Changed statement tags in CONTINUE-statement
    98 ! 
     101!
    99102! 1367 2014-04-23 15:18:30Z witha
    100103! Bugfix: module kinds must be used
     
    103106! +alloc_factor, + min_nr_particle
    104107! -dt_sort_particles, -maximum_number_of_particles
    105 ! 
     108!
    106109! 1340 2014-03-25 19:45:13Z kanani
    107110! REAL constants defined as wp-kinds
     
    138141!> software packages which are used optionally in the run.
    139142!>
    140 !> @todo Perform all actions in the respective submodules and remove 
     143!> @todo Perform all actions in the respective submodules and remove
    141144!>       package_parin
    142145!------------------------------------------------------------------------------!
    143146 SUBROUTINE package_parin
    144  
     147
    145148
    146149    USE control_parameters,                                                    &
     
    165168
    166169    USE particle_attributes,                                                   &
    167         ONLY:  alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,         &
    168                collision_algorithm, collision_kernel,                          &
     170        ONLY:  aero_type, aero_weight, alloc_factor, bc_par_b, bc_par_lr,      &
     171               bc_par_ns, bc_par_t, collision_kernel,                          &
    169172               curvature_solution_effects, deallocate_memory, density_ratio,   &
    170173               dissipation_classes, dt_min_part, dt_prel,                      &
    171174               dt_write_particle_data, end_time_prel, initial_weighting_factor,&
    172                init_aerosol_probabilistic, max_number_particles_per_gridbox,   &
    173                merging, min_nr_particle, monodisperse_aerosols, n1, n2, n3,    &
     175               log_sigma, max_number_particles_per_gridbox,                    &
     176               merging, min_nr_particle, na,                                   &
    174177               number_concentration, number_of_particle_groups,                &
    175178               number_particles_per_gridbox, particles_per_point,              &
     
    177180               psb, psl, psn, psr, pss, pst, radius, radius_classes,           &
    178181               radius_merge, radius_split, random_start_position,              &
    179                read_particles_from_restartfile, rm1, rm2, rm3,                 &
    180                seed_follows_topography, splitting, splitting_factor,           & 
     182               read_particles_from_restartfile, rm,                            &
     183               seed_follows_topography, splitting, splitting_factor,           &
    181184               splitting_factor_max, splitting_function, splitting_mode,       &
    182                step_dealloc, s1, s2, s3, use_sgs_for_particles,                &
     185               step_dealloc, use_sgs_for_particles,                            &
    183186               vertical_particle_advection, weight_factor_merge,               &
    184                weight_factor_split, write_particle_statistics               
     187               weight_factor_split, write_particle_statistics
    185188
    186189    IMPLICIT NONE
     
    206209                                  vc_size_y, vc_size_z
    207210
    208     NAMELIST /particles_par/      alloc_factor, bc_par_b, bc_par_lr,           &
    209                                   bc_par_ns, bc_par_t, collision_algorithm,    &
     211    NAMELIST /particles_par/      aero_type, aero_weight, alloc_factor,        &
     212                                  bc_par_b, bc_par_lr,                         &
     213                                  bc_par_ns, bc_par_t,                         &
    210214                                  collision_kernel, curvature_solution_effects,&
    211215                                  deallocate_memory, density_ratio,            &
     
    214218                                  dt_write_particle_data,                      &
    215219                                  end_time_prel, initial_weighting_factor,     &
    216                                   init_aerosol_probabilistic,                  &
     220                                  log_sigma,                                   &
    217221                                  max_number_particles_per_gridbox, merging,   &
    218                                   min_nr_particle, monodisperse_aerosols,      &
    219                                   n1, n2, n3, number_concentration,            &
     222                                  min_nr_particle,                             &
     223                                  na, number_concentration,                    &
    220224                                  number_of_particle_groups,                   &
    221225                                  number_particles_per_gridbox,                &
     
    224228                                  particle_maximum_age, pdx, pdy, pdz, psb,    &
    225229                                  psl, psn, psr, pss, pst, radius,             &
    226                                   radius_classes, radius_merge, radius_split,  &   
     230                                  radius_classes, radius_merge, radius_split,  &
    227231                                  random_start_position,                       &
    228                                   read_particles_from_restartfile,             &
    229                                   rm1, rm2, rm3,                               &
     232                                  read_particles_from_restartfile, rm,         &
    230233                                  seed_follows_topography,                     &
    231234                                  splitting, splitting_factor,                 &
    232235                                  splitting_factor_max, splitting_function,    &
    233                                   splitting_mode, step_dealloc, s1, s2, s3,    &
     236                                  splitting_mode, step_dealloc,                &
    234237                                  use_sgs_for_particles,                       &
    235238                                  vertical_particle_advection,                 &
Note: See TracChangeset for help on using the changeset viewer.