source: palm/trunk/SOURCE/salsa_mod.f90 @ 4109

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

  • Property svn:keywords set to Id
File size: 505.3 KB
Line 
1!> @file salsa_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM-4U.
4!
5! PALM-4U is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2018-2019 University of Helsinki
18! Copyright 1997-2019 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23! Pass integer flag array as well as boundary flags to WS scalar advection
24! routine
25!
26! Former revisions:
27! -----------------
28! $Id: salsa_mod.f90 4109 2019-07-22 17:00:34Z suehring $
29! Slightly revise setting of boundary conditions at horizontal walls, use
30! data-structure offset index instead of pre-calculate it for each facing
31!
32! 4079 2019-07-09 18:04:41Z suehring
33! Application of monotonic flux limiter for the vertical scalar advection
34! up to the topography top (only for the cache-optimized version at the
35! moment).
36!
37! 4069 2019-07-01 14:05:51Z Giersch
38! Masked output running index mid has been introduced as a local variable to
39! avoid runtime error (Loop variable has been modified) in time_integration
40!
41! 4058 2019-06-27 15:25:42Z knoop
42! Bugfix: to_be_resorted was uninitialized in case of s_H2O in 3d_data_averaging
43!
44! 4012 2019-05-31 15:19:05Z monakurppa
45! Merge salsa branch to trunk. List of changes:
46! - Error corrected in distr_update that resulted in the aerosol number size
47!   distribution not converging if the concentration was nclim.
48! - Added a separate output for aerosol liquid water (s_H2O)
49! - aerosol processes for a size bin are now calculated only if the aerosol
50!   number of concentration of that bin is > 2*nclim
51! - An initialisation error in the subroutine "deposition" corrected and the
52!   subroutine reformatted.
53! - stuff from salsa_util_mod.f90 moved into salsa_mod.f90
54! - calls for closing the netcdf input files added
55!
56! 3956 2019-05-07 12:32:52Z monakurppa
57! - Conceptual bug in depo_surf correct for urban and land surface model
58! - Subroutine salsa_tendency_ij optimized.
59! - Interfaces salsa_non_advective_processes and salsa_exchange_horiz_bounds
60!   created. These are now called in module_interface.
61!   salsa_exchange_horiz_bounds after calling salsa_driver only when needed
62!   (i.e. every dt_salsa).
63!
64! 3924 2019-04-23 09:33:06Z monakurppa
65! Correct a bug introduced by the previous update.
66!
67! 3899 2019-04-16 14:05:27Z monakurppa
68! - remove unnecessary error / location messages
69! - corrected some error message numbers
70! - allocate source arrays only if emissions or dry deposition is applied.
71!
72! 3885 2019-04-11 11:29:34Z kanani
73! Changes related to global restructuring of location messages and introduction
74! of additional debug messages
75!
76! 3876 2019-04-08 18:41:49Z knoop
77! Introduced salsa_actions module interface
78!
79! 3871 2019-04-08 14:38:39Z knoop
80! Major changes in formatting, performance and data input structure (see branch
81! the history for details)
82! - Time-dependent emissions enabled: lod=1 for yearly PM emissions that are
83!   normalised depending on the time, and lod=2 for preprocessed emissions
84!   (similar to the chemistry module).
85! - Additionally, 'uniform' emissions allowed. This emission is set constant on
86!   all horisontal upward facing surfaces and it is created based on parameters
87!   surface_aerosol_flux, aerosol_flux_dpg/sigmag/mass_fracs_a/mass_fracs_b.
88! - All emissions are now implemented as surface fluxes! No 3D sources anymore.
89! - Update the emission information by calling salsa_emission_update if
90!   skip_time_do_salsa >= time_since_reference_point and
91!   next_aero_emission_update <= time_since_reference_point
92! - Aerosol background concentrations read from PIDS_DYNAMIC. The vertical grid
93!   must match the one applied in the model.
94! - Gas emissions and background concentrations can be also read in in salsa_mod
95!   if the chemistry module is not applied.
96! - In deposition, information on the land use type can be now imported from
97!   the land use model
98! - Use SI units in PARIN, i.e. n_lognorm given in #/m3 and dpg in metres.
99! - Apply 100 character line limit
100! - Change all variable names from capital to lowercase letter
101! - Change real exponents to integer if possible. If not, precalculate the value
102!   value of exponent
103! - Rename in1a to start_subrange_1a, fn2a to end_subrange_1a etc.
104! - Rename nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and ngast -->
105!   ngases_salsa
106! - Rename ibc to index_bc, idu to index_du etc.
107! - Renamed loop indices b, c and sg to ib, ic and ig
108! - run_salsa subroutine removed
109! - Corrected a bud in salsa_driver: falsely applied ino instead of inh
110! - Call salsa_tendency within salsa_prognostic_equations which is called in
111!   module_interface_mod instead of prognostic_equations_mod
112! - Removed tailing white spaces and unused variables
113! - Change error message to start by PA instead of SA
114!
115! 3833 2019-03-28 15:04:04Z forkel
116! added USE chem_gasphase_mod for nvar, nspec and spc_names
117!
118! 3787 2019-03-07 08:43:54Z raasch
119! unused variables removed
120!
121! 3780 2019-03-05 11:19:45Z forkel
122! unused variable for file index removed from rrd-subroutines parameter list
123!
124! 3685 2019-01-21 01:02:11Z knoop
125! Some interface calls moved to module_interface + cleanup
126!
127! 3655 2019-01-07 16:51:22Z knoop
128! Implementation of the PALM module interface
129!
130! 3636 2018-12-19 13:48:34Z raasch
131! nopointer option removed
132!
133! 3630 2018-12-17 11:04:17Z knoop
134! - Moved the control parameter "salsa" from salsa_mod.f90 to control_parameters
135! - Updated salsa_rrd_local and salsa_wrd_local
136! - Add target attribute
137! - Revise initialization in case of restarts
138! - Revise masked data output
139!
140! 3582 2018-11-29 19:16:36Z suehring
141! missing comma separator inserted
142!
143! 3483 2018-11-02 14:19:26Z raasch
144! bugfix: directives added to allow compilation without netCDF
145!
146! 3481 2018-11-02 09:14:13Z raasch
147! temporary variable cc introduced to circumvent a possible Intel18 compiler bug
148! related to contiguous/non-contguous pointer/target attributes
149!
150! 3473 2018-10-30 20:50:15Z suehring
151! NetCDF input routine renamed
152!
153! 3467 2018-10-30 19:05:21Z suehring
154! Initial revision
155!
156! 3412 2018-10-24 07:25:57Z monakurppa
157!
158! Authors:
159! --------
160! @author Mona Kurppa (University of Helsinki)
161!
162!
163! Description:
164! ------------
165!> Sectional aerosol module for large scale applications SALSA
166!> (Kokkola et al., 2008, ACP 8, 2469-2483). Solves the aerosol number and mass
167!> concentration as well as chemical composition. Includes aerosol dynamic
168!> processes: nucleation, condensation/evaporation of vapours, coagulation and
169!> deposition on tree leaves, ground and roofs.
170!> Implementation is based on formulations implemented in UCLALES-SALSA except
171!> for deposition which is based on parametrisations by Zhang et al. (2001,
172!> Atmos. Environ. 35, 549-560) or Petroff&Zhang (2010, Geosci. Model Dev. 3,
173!> 753-769)
174!>
175!> @todo Apply information from emission_stack_height to lift emission sources
176!> @todo emission mode "parameterized", i.e. based on street type
177!> @todo Allow insoluble emissions
178!> @todo Apply flux limiter in prognostic equations
179!------------------------------------------------------------------------------!
180 MODULE salsa_mod
181
182    USE basic_constants_and_equations_mod,                                                         &
183        ONLY:  c_p, g, p_0, pi, r_d
184
185    USE chem_gasphase_mod,                                                                         &
186        ONLY:  nspec, nvar, spc_names
187
188    USE chem_modules,                                                                              &
189        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, chem_species
190
191    USE control_parameters
192
193    USE indices,                                                                                   &
194        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nz, nzt, wall_flags_0
195
196    USE kinds
197
198    USE netcdf_data_input_mod,                                                                     &
199        ONLY:  chem_emis_att_type, chem_emis_val_type
200
201    USE pegrid
202
203    USE statistics,                                                                                &
204        ONLY:  sums_salsa_ws_l
205
206    IMPLICIT NONE
207!
208!-- SALSA constants:
209!
210!-- Local constants:
211    INTEGER(iwp), PARAMETER ::  luc_urban = 15     !< default landuse type for urban
212    INTEGER(iwp), PARAMETER ::  ngases_salsa  = 5  !< total number of gaseous tracers:
213                                                   !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV
214                                                   !< (non-volatile OC), 5 = OCSV (semi-volatile)
215    INTEGER(iwp), PARAMETER ::  nmod = 7     !< number of modes for initialising the aerosol size
216                                             !< distribution
217    INTEGER(iwp), PARAMETER ::  nreg = 2     !< Number of main size subranges
218    INTEGER(iwp), PARAMETER ::  maxspec = 7  !< Max. number of aerosol species
219    INTEGER(iwp), PARAMETER ::  season = 1   !< For dry depostion by Zhang et al.: 1 = summer,
220                                             !< 2 = autumn (no harvest yet), 3 = late autumn
221                                             !< (already frost), 4 = winter, 5 = transitional spring
222!
223!-- Universal constants
224    REAL(wp), PARAMETER ::  abo    = 1.380662E-23_wp   !< Boltzmann constant (J/K)
225    REAL(wp), PARAMETER ::  alv    = 2.260E+6_wp       !< latent heat for H2O
226                                                       !< vaporisation (J/kg)
227    REAL(wp), PARAMETER ::  alv_d_rv  = 4896.96865_wp  !< alv / rv
228    REAL(wp), PARAMETER ::  am_airmol = 4.8096E-26_wp  !< Average mass of one air
229                                                       !< molecule (Jacobson,
230                                                       !< 2005, Eq. 2.3)
231    REAL(wp), PARAMETER ::  api6   = 0.5235988_wp      !< pi / 6
232    REAL(wp), PARAMETER ::  argas  = 8.314409_wp       !< Gas constant (J/(mol K))
233    REAL(wp), PARAMETER ::  argas_d_cpd = 8.281283865E-3_wp  !< argas per cpd
234    REAL(wp), PARAMETER ::  avo    = 6.02214E+23_wp    !< Avogadro constant (1/mol)
235    REAL(wp), PARAMETER ::  d_sa   = 5.539376964394570E-10_wp  !< diameter of condensing sulphuric
236                                                               !< acid molecule (m)
237    REAL(wp), PARAMETER ::  for_ppm_to_nconc =  7.243016311E+16_wp !< ppm * avo / R (K/(Pa*m3))
238    REAL(wp), PARAMETER ::  epsoc  = 0.15_wp          !< water uptake of organic
239                                                      !< material
240    REAL(wp), PARAMETER ::  mclim  = 1.0E-23_wp       !< mass concentration min limit (kg/m3)
241    REAL(wp), PARAMETER ::  n3     = 158.79_wp        !< Number of H2SO4 molecules in 3 nm cluster
242                                                      !< if d_sa=5.54e-10m
243    REAL(wp), PARAMETER ::  nclim  = 1.0_wp           !< number concentration min limit (#/m3)
244    REAL(wp), PARAMETER ::  surfw0 = 0.073_wp         !< surface tension of water at 293 K (J/m2)
245!
246!-- Molar masses in kg/mol
247    REAL(wp), PARAMETER ::  ambc     = 12.0E-3_wp     !< black carbon (BC)
248    REAL(wp), PARAMETER ::  amdair   = 28.970E-3_wp   !< dry air
249    REAL(wp), PARAMETER ::  amdu     = 100.E-3_wp     !< mineral dust
250    REAL(wp), PARAMETER ::  amh2o    = 18.0154E-3_wp  !< H2O
251    REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp    !< H2SO4
252    REAL(wp), PARAMETER ::  amhno3   = 63.01E-3_wp    !< HNO3
253    REAL(wp), PARAMETER ::  amn2o    = 44.013E-3_wp   !< N2O
254    REAL(wp), PARAMETER ::  amnh3    = 17.031E-3_wp   !< NH3
255    REAL(wp), PARAMETER ::  amo2     = 31.9988E-3_wp  !< O2
256    REAL(wp), PARAMETER ::  amo3     = 47.998E-3_wp   !< O3
257    REAL(wp), PARAMETER ::  amoc     = 150.E-3_wp     !< organic carbon (OC)
258    REAL(wp), PARAMETER ::  amss     = 58.44E-3_wp    !< sea salt (NaCl)
259!
260!-- Densities in kg/m3
261    REAL(wp), PARAMETER ::  arhobc     = 2000.0_wp  !< black carbon
262    REAL(wp), PARAMETER ::  arhodu     = 2650.0_wp  !< mineral dust
263    REAL(wp), PARAMETER ::  arhoh2o    = 1000.0_wp  !< H2O
264    REAL(wp), PARAMETER ::  arhoh2so4  = 1830.0_wp  !< SO4
265    REAL(wp), PARAMETER ::  arhohno3   = 1479.0_wp  !< HNO3
266    REAL(wp), PARAMETER ::  arhonh3    = 1530.0_wp  !< NH3
267    REAL(wp), PARAMETER ::  arhooc     = 2000.0_wp  !< organic carbon
268    REAL(wp), PARAMETER ::  arhoss     = 2165.0_wp  !< sea salt (NaCl)
269!
270!-- Volume of molecule in m3/#
271    REAL(wp), PARAMETER ::  amvh2o   = amh2o /avo / arhoh2o      !< H2O
272    REAL(wp), PARAMETER ::  amvh2so4 = amh2so4 / avo / arhoh2so4 !< SO4
273    REAL(wp), PARAMETER ::  amvhno3  = amhno3 / avo / arhohno3   !< HNO3
274    REAL(wp), PARAMETER ::  amvnh3   = amnh3 / avo / arhonh3     !< NH3
275    REAL(wp), PARAMETER ::  amvoc    = amoc / avo / arhooc       !< OC
276    REAL(wp), PARAMETER ::  amvss    = amss / avo / arhoss       !< sea salt
277!
278!-- Constants for the dry deposition model by Petroff and Zhang (2010):
279!-- obstacle characteristic dimension "L" (cm) (plane obstacle by default) and empirical constants
280!-- C_B, C_IN, C_IM, beta_IM and C_IT for each land use category (15, as in Zhang et al. (2001))
281    REAL(wp), DIMENSION(1:15), PARAMETER :: l_p10 = &
282        (/0.15,  4.0,   0.15,  3.0,   3.0,   0.5,   3.0,   -99., 0.5,  2.0,  1.0,   -99., -99., -99., 3.0/)
283    REAL(wp), DIMENSION(1:15), PARAMETER :: c_b_p10 = &
284        (/0.887, 1.262, 0.887, 1.262, 1.262, 0.996, 0.996, -99., 0.7,  0.93, 0.996, -99., -99., -99., 1.262/)
285    REAL(wp), DIMENSION(1:15), PARAMETER :: c_in_p10 = &
286        (/0.81,  0.216, 0.81,  0.216, 0.216, 0.191, 0.162, -99., 0.7,  0.14, 0.162, -99., -99., -99., 0.216/)
287    REAL(wp), DIMENSION(1:15), PARAMETER :: c_im_p10 = &
288        (/0.162, 0.13,  0.162, 0.13,  0.13,  0.191, 0.081, -99., 0.191,0.086,0.081, -99., -99., -99., 0.13/)
289    REAL(wp), DIMENSION(1:15), PARAMETER :: beta_im_p10 = &
290        (/0.6,   0.47,  0.6,   0.47,  0.47,  0.47,  0.47,  -99., 0.6,  0.47, 0.47,  -99., -99., -99., 0.47/)
291    REAL(wp), DIMENSION(1:15), PARAMETER :: c_it_p10 = &
292        (/0.0,   0.056, 0.0,   0.056, 0.056, 0.042, 0.056, -99., 0.042,0.014,0.056, -99., -99., -99., 0.056/)
293!
294!-- Constants for the dry deposition model by Zhang et al. (2001):
295!-- empirical constants "alpha" and "gamma" and characteristic radius "A" for
296!-- each land use category (15) and season (5)
297    REAL(wp), DIMENSION(1:15), PARAMETER :: alpha_z01 = &
298        (/1.0, 0.6, 1.1, 0.8, 0.8, 1.2, 1.2, 50.0, 50.0, 1.3, 2.0, 50.0, 100.0, 100.0, 1.5/)
299    REAL(wp), DIMENSION(1:15), PARAMETER :: gamma_z01 = &
300        (/0.56, 0.58, 0.56, 0.56, 0.56, 0.54, 0.54, 0.54, 0.54, 0.54, 0.54, 0.54, 0.50, 0.50, 0.56/)
301    REAL(wp), DIMENSION(1:15,1:5), PARAMETER :: A_z01 =  RESHAPE( (/& 
302         2.0, 5.0, 2.0,  5.0, 5.0, 2.0, 2.0, -99., -99., 10.0, 10.0, -99., -99., -99., 10.0,&  ! SC1
303         2.0, 5.0, 2.0,  5.0, 5.0, 2.0, 2.0, -99., -99., 10.0, 10.0, -99., -99., -99., 10.0,&  ! SC2
304         2.0, 5.0, 5.0, 10.0, 5.0, 5.0, 5.0, -99., -99., 10.0, 10.0, -99., -99., -99., 10.0,&  ! SC3
305         2.0, 5.0, 5.0, 10.0, 5.0, 5.0, 5.0, -99., -99., 10.0, 10.0, -99., -99., -99., 10.0,&  ! SC4
306         2.0, 5.0, 2.0,  5.0, 5.0, 2.0, 2.0, -99., -99., 10.0, 10.0, -99., -99., -99., 10.0 &  ! SC5
307                                                           /), (/ 15, 5 /) )
308!-- Land use categories (based on Z01 but the same applies here also for P10):
309!-- 1 = evergreen needleleaf trees,
310!-- 2 = evergreen broadleaf trees,
311!-- 3 = deciduous needleleaf trees,
312!-- 4 = deciduous broadleaf trees,
313!-- 5 = mixed broadleaf and needleleaf trees (deciduous broadleaf trees for P10),
314!-- 6 = grass (short grass for P10),
315!-- 7 = crops, mixed farming,
316!-- 8 = desert,
317!-- 9 = tundra,
318!-- 10 = shrubs and interrupted woodlands (thorn shrubs for P10),
319!-- 11 = wetland with plants (long grass for P10)
320!-- 12 = ice cap and glacier,
321!-- 13 = inland water (inland lake for P10)
322!-- 14 = ocean (water for P10),
323!-- 15 = urban
324!
325!-- SALSA variables:
326    CHARACTER(LEN=20)  ::  bc_salsa_b = 'neumann'                 !< bottom boundary condition
327    CHARACTER(LEN=20)  ::  bc_salsa_t = 'neumann'                 !< top boundary condition
328    CHARACTER(LEN=20)  ::  depo_pcm_par = 'zhang2001'             !< or 'petroff2010'
329    CHARACTER(LEN=20)  ::  depo_pcm_type = 'deciduous_broadleaf'  !< leaf type
330    CHARACTER(LEN=20)  ::  depo_surf_par = 'zhang2001'            !< or 'petroff2010'
331    CHARACTER(LEN=100) ::  input_file_dynamic = 'PIDS_DYNAMIC'    !< file name for dynamic input
332    CHARACTER(LEN=100) ::  input_file_salsa   = 'PIDS_SALSA'      !< file name for emission data
333    CHARACTER(LEN=20)  ::  salsa_emission_mode = 'no_emission'    !< 'no_emission', 'uniform',
334                                                                  !< 'parameterized', 'read_from_file'
335
336    CHARACTER(LEN=20), DIMENSION(4) ::  decycle_method =                                           &
337                                                 (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
338                                     !< Decycling method at horizontal boundaries
339                                     !< 1=left, 2=right, 3=south, 4=north
340                                     !< dirichlet = initial profiles for the ghost and first 3 layers
341                                     !< neumann = zero gradient
342
343    CHARACTER(LEN=3), DIMENSION(maxspec) ::  listspec = &  !< Active aerosols
344                                   (/'SO4','   ','   ','   ','   ','   ','   '/)
345
346    INTEGER(iwp) ::  depo_pcm_par_num = 1   !< parametrisation type: 1=zhang2001, 2=petroff2010
347    INTEGER(iwp) ::  depo_pcm_type_num = 0  !< index for the dry deposition type on the plant canopy
348    INTEGER(iwp) ::  depo_surf_par_num = 1  !< parametrisation type: 1=zhang2001, 2=petroff2010
349    INTEGER(iwp) ::  dots_salsa = 0         !< starting index for salsa-timeseries
350    INTEGER(iwp) ::  end_subrange_1a = 1    !< last index for bin subrange 1a
351    INTEGER(iwp) ::  end_subrange_2a = 1    !< last index for bin subrange 2a
352    INTEGER(iwp) ::  end_subrange_2b = 1    !< last index for bin subrange 2b
353    INTEGER(iwp) ::  ibc_salsa_b            !< index for the bottom boundary condition
354    INTEGER(iwp) ::  ibc_salsa_t            !< index for the top boundary condition
355    INTEGER(iwp) ::  index_bc  = -1         !< index for black carbon (BC)
356    INTEGER(iwp) ::  index_du  = -1         !< index for dust
357    INTEGER(iwp) ::  index_nh  = -1         !< index for NH3
358    INTEGER(iwp) ::  index_no  = -1         !< index for HNO3
359    INTEGER(iwp) ::  index_oc  = -1         !< index for organic carbon (OC)
360    INTEGER(iwp) ::  index_so4 = -1         !< index for SO4 or H2SO4
361    INTEGER(iwp) ::  index_ss  = -1         !< index for sea salt
362    INTEGER(iwp) ::  init_aerosol_type = 0  !< Initial size distribution type
363                                            !< 0 = uniform (read from PARIN)
364                                            !< 1 = read vertical profile of the mode number
365                                            !<     concentration from an input file
366    INTEGER(iwp) ::  init_gases_type = 0    !< Initial gas concentration type
367                                            !< 0 = uniform (read from PARIN)
368                                            !< 1 = read vertical profile from an input file
369    INTEGER(iwp) ::  lod_gas_emissions = 0  !< level of detail of the gaseous emission data
370    INTEGER(iwp) ::  nbins_aerosol = 1      !< total number of size bins
371    INTEGER(iwp) ::  ncc   = 1              !< number of chemical components used
372    INTEGER(iwp) ::  ncomponents_mass = 1   !< total number of chemical compounds (ncc+1)
373                                            !< if particle water is advected)
374    INTEGER(iwp) ::  nj3 = 1                !< J3 parametrization (nucleation)
375                                            !< 1 = condensational sink (Kerminen&Kulmala, 2002)
376                                            !< 2 = coagulational sink (Lehtinen et al. 2007)
377                                            !< 3 = coagS+self-coagulation (Anttila et al. 2010)
378    INTEGER(iwp) ::  nsnucl = 0             !< Choice of the nucleation scheme:
379                                            !< 0 = off
380                                            !< 1 = binary nucleation
381                                            !< 2 = activation type nucleation
382                                            !< 3 = kinetic nucleation
383                                            !< 4 = ternary nucleation
384                                            !< 5 = nucleation with ORGANICs
385                                            !< 6 = activation type of nucleation with H2SO4+ORG
386                                            !< 7 = heteromolecular nucleation with H2SO4*ORG
387                                            !< 8 = homomolecular nucleation of H2SO4
388                                            !<     + heteromolecular nucleation with H2SO4*ORG
389                                            !< 9 = homomolecular nucleation of H2SO4 and ORG
390                                            !<     + heteromolecular nucleation with H2SO4*ORG
391    INTEGER(iwp) ::  start_subrange_1a = 1  !< start index for bin subranges: subrange 1a
392    INTEGER(iwp) ::  start_subrange_2a = 1  !<                                subrange 2a
393    INTEGER(iwp) ::  start_subrange_2b = 1  !<                                subrange 2b
394
395    INTEGER(iwp), DIMENSION(nreg) ::  nbin = (/ 3, 7/)  !< Number of size bins per subrange: 1 & 2
396
397    INTEGER(iwp), DIMENSION(ngases_salsa) ::  gas_index_chem = &
398                                                 (/ 1, 1, 1, 1, 1/)  !< gas indices in chemistry_model_mod
399                                                 !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV, 5 = OCSV
400    INTEGER(iwp), DIMENSION(ngases_salsa) ::  emission_index_chem  !< gas indices in the gas emission file
401
402    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  k_topo_top  !< vertical index of the topography top
403!
404!-- SALSA switches:
405    LOGICAL ::  advect_particle_water   = .TRUE.   !< Advect water concentration of particles
406    LOGICAL ::  decycle_lr              = .FALSE.  !< Undo cyclic boundaries: left and right
407    LOGICAL ::  decycle_ns              = .FALSE.  !< Undo cyclic boundaries: north and south
408    LOGICAL ::  include_emission        = .FALSE.  !< Include or not emissions
409    LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
410    LOGICAL ::  nest_salsa              = .FALSE.  !< Apply nesting for salsa
411    LOGICAL ::  no_insoluble            = .FALSE.  !< Exclude insoluble chemical components
412    LOGICAL ::  read_restart_data_salsa = .FALSE.  !< Read restart data for salsa
413    LOGICAL ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA from
414                                                   !< the chemistry model
415    LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Enhancement of coagulation kernel by van der
416                                                   !< Waals and viscous forces
417    LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
418!
419!-- Process switches: nl* is read from the NAMELIST and is NOT changed.
420!--                   ls* is the switch used and will get the value of nl*
421!--                       except for special circumstances (spinup period etc.)
422    LOGICAL ::  nlcoag       = .FALSE.  !< Coagulation master switch
423    LOGICAL ::  lscoag       = .FALSE.  !<
424    LOGICAL ::  nlcnd        = .FALSE.  !< Condensation master switch
425    LOGICAL ::  lscnd        = .FALSE.  !<
426    LOGICAL ::  nlcndgas     = .FALSE.  !< Condensation of precursor gases
427    LOGICAL ::  lscndgas     = .FALSE.  !<
428    LOGICAL ::  nlcndh2oae   = .FALSE.  !< Condensation of H2O on aerosol
429    LOGICAL ::  lscndh2oae   = .FALSE.  !< particles (FALSE -> equilibrium calc.)
430    LOGICAL ::  nldepo       = .FALSE.  !< Deposition master switch
431    LOGICAL ::  lsdepo       = .FALSE.  !<
432    LOGICAL ::  nldepo_surf  = .FALSE.  !< Deposition on vegetation master switch
433    LOGICAL ::  lsdepo_surf  = .FALSE.  !<
434    LOGICAL ::  nldepo_pcm   = .FALSE.  !< Deposition on walls master switch
435    LOGICAL ::  lsdepo_pcm   = .FALSE.  !<
436    LOGICAL ::  nldistupdate = .TRUE.   !< Size distribution update master switch
437    LOGICAL ::  lsdistupdate = .FALSE.  !<
438    LOGICAL ::  lspartition  = .FALSE.  !< Partition of HNO3 and NH3
439
440    REAL(wp) ::  act_coeff = 1.0E-7_wp               !< Activation coefficient
441    REAL(wp) ::  dt_salsa  = 0.00001_wp              !< Time step of SALSA
442    REAL(wp) ::  h2so4_init = nclim                  !< Init value for sulphuric acid gas
443    REAL(wp) ::  hno3_init  = nclim                  !< Init value for nitric acid gas
444    REAL(wp) ::  last_salsa_time = 0.0_wp            !< previous salsa call
445    REAL(wp) ::  next_aero_emission_update = 0.0_wp  !< previous emission update
446    REAL(wp) ::  next_gas_emission_update = 0.0_wp   !< previous emission update
447    REAL(wp) ::  nf2a = 1.0_wp                       !< Number fraction allocated to 2a-bins
448    REAL(wp) ::  nh3_init  = nclim                   !< Init value for ammonia gas
449    REAL(wp) ::  ocnv_init = nclim                   !< Init value for non-volatile organic gases
450    REAL(wp) ::  ocsv_init = nclim                   !< Init value for semi-volatile organic gases
451    REAL(wp) ::  rhlim = 1.20_wp                     !< RH limit in %/100. Prevents unrealistical RH
452    REAL(wp) ::  skip_time_do_salsa = 0.0_wp         !< Starting time of SALSA (s)
453!
454!-- Initial log-normal size distribution: mode diameter (dpg, metres),
455!-- standard deviation (sigmag) and concentration (n_lognorm, #/m3)
456    REAL(wp), DIMENSION(nmod) ::  dpg   = &
457                     (/1.3E-8_wp, 5.4E-8_wp, 8.6E-7_wp, 2.0E-7_wp, 2.0E-7_wp, 2.0E-7_wp, 2.0E-7_wp/)
458    REAL(wp), DIMENSION(nmod) ::  sigmag  = &
459                                        (/1.8_wp, 2.16_wp, 2.21_wp, 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/)
460    REAL(wp), DIMENSION(nmod) ::  n_lognorm = &
461                             (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
462!
463!-- Initial mass fractions / chemical composition of the size distribution
464    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_a = & !< mass fractions between
465             (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for A bins
466    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_b = & !< mass fractions between
467             (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for B bins
468    REAL(wp), DIMENSION(nreg+1) ::  reglim = & !< Min&max diameters of size subranges
469                                 (/ 3.0E-9_wp, 5.0E-8_wp, 1.0E-5_wp/)
470!
471!-- Initial log-normal size distribution: mode diameter (dpg, metres), standard deviation (sigmag)
472!-- concentration (n_lognorm, #/m3) and mass fractions of all chemical components (listed in
473!-- listspec) for both a (soluble) and b (insoluble) bins.
474    REAL(wp), DIMENSION(nmod) ::  aerosol_flux_dpg   = &
475                     (/1.3E-8_wp, 5.4E-8_wp, 8.6E-7_wp, 2.0E-7_wp, 2.0E-7_wp, 2.0E-7_wp, 2.0E-7_wp/)
476    REAL(wp), DIMENSION(nmod) ::  aerosol_flux_sigmag  = &
477                                        (/1.8_wp, 2.16_wp, 2.21_wp, 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/)
478    REAL(wp), DIMENSION(maxspec) ::  aerosol_flux_mass_fracs_a = &
479                                                               (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
480    REAL(wp), DIMENSION(maxspec) ::  aerosol_flux_mass_fracs_b = &
481                                                               (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
482    REAL(wp), DIMENSION(nmod) ::  surface_aerosol_flux = &
483                                 (/1.0E+8_wp, 1.0E+9_wp, 1.0E+5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
484
485    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bin_low_limits     !< to deliver information about
486                                                               !< the lower diameters per bin
487    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bc_am_t_val        !< vertical gradient of: aerosol mass
488    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bc_an_t_val        !< of: aerosol number
489    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bc_gt_t_val        !< salsa gases near domain top
490    REAL(wp), DIMENSION(:), ALLOCATABLE ::  gas_emission_time  !< Time array in gas emission data (s)
491    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nsect              !< Background number concentrations
492    REAL(wp), DIMENSION(:), ALLOCATABLE ::  massacc            !< Mass accomodation coefficients
493!
494!-- SALSA derived datatypes:
495!
496!-- Component index
497    TYPE component_index
498       CHARACTER(len=3), ALLOCATABLE ::  comp(:)  !< Component name
499       INTEGER(iwp) ::  ncomp  !< Number of components
500       INTEGER(iwp), ALLOCATABLE ::  ind(:)  !< Component index
501    END TYPE component_index
502!
503!-- For matching LSM and USM surface types and the deposition module surface types
504    TYPE match_surface
505       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_lupg  !< index for pavement / green roofs
506       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luvw  !< index for vegetation / walls
507       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luww  !< index for water / windows
508    END TYPE match_surface
509!
510!-- Aerosol emission data attributes
511    TYPE salsa_emission_attribute_type
512
513       CHARACTER(LEN=25) ::   units
514
515       CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE ::   cat_name    !<
516       CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE ::   cc_name     !<
517       CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE ::   unit_time   !<
518       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names   !<
519
520       INTEGER(iwp) ::  lod = 0            !< level of detail
521       INTEGER(iwp) ::  nbins = 10         !< number of aerosol size bins
522       INTEGER(iwp) ::  ncat  = 0          !< number of emission categories
523       INTEGER(iwp) ::  ncc   = 7          !< number of aerosol chemical components
524       INTEGER(iwp) ::  nhoursyear = 0     !< number of hours: HOURLY mode
525       INTEGER(iwp) ::  nmonthdayhour = 0  !< number of month days and hours: MDH mode
526       INTEGER(iwp) ::  num_vars           !< number of variables
527       INTEGER(iwp) ::  nt  = 0            !< number of time steps
528       INTEGER(iwp) ::  nz  = 0            !< number of vertical levels
529       INTEGER(iwp) ::  tind               !< time index for reference time in salsa emission data
530
531       INTEGER(iwp), DIMENSION(maxspec) ::  cc_input_to_model   !<
532
533       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  cat_index  !< Index of emission categories
534       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  cc_index   !< Index of chemical components
535
536       REAL(wp) ::  conversion_factor  !< unit conversion factor for aerosol emissions
537
538       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dmid         !< mean diameters of size bins (m)
539       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rho          !< average density (kg/m3)
540       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time         !< time (s)
541       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_factor  !< emission time factor
542       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z            !< height (m)
543
544       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  etf  !< emission time factor
545       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: stack_height
546
547    END TYPE salsa_emission_attribute_type
548!
549!-- The default size distribution and mass composition per emission category:
550!-- 1 = traffic, 2 = road dust, 3 = wood combustion, 4 = other
551!-- Mass fractions: H2SO4, OC, BC, DU, SS, HNO3, NH3
552    TYPE salsa_emission_mode_type
553
554       INTEGER(iwp) ::  ndm = 3  !< number of default modes
555       INTEGER(iwp) ::  ndc = 4  !< number of default categories
556
557       CHARACTER(LEN=25), DIMENSION(1:4) ::  cat_name_table = (/'traffic exhaust', &
558                                                                'road dust      ', &
559                                                                'wood combustion', &
560                                                                'other          '/)
561
562       INTEGER(iwp), DIMENSION(1:4) ::  cat_input_to_model   !<
563
564       REAL(wp), DIMENSION(1:3) ::  dpg_table = (/ 13.5E-9_wp, 1.4E-6_wp, 5.4E-8_wp/)  !<
565       REAL(wp), DIMENSION(1:3) ::  ntot_table  !<
566       REAL(wp), DIMENSION(1:3) ::  sigmag_table = (/ 1.6_wp, 1.4_wp, 1.7_wp /)  !<
567
568       REAL(wp), DIMENSION(1:maxspec,1:4) ::  mass_frac_table = &  !<
569          RESHAPE( (/ 0.04_wp, 0.48_wp, 0.48_wp, 0.0_wp,  0.0_wp, 0.0_wp, 0.0_wp, &
570                      0.0_wp,  0.05_wp, 0.0_wp,  0.95_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
571                      0.0_wp,  0.5_wp,  0.5_wp,  0.0_wp,  0.0_wp, 0.0_wp, 0.0_wp, &
572                      0.0_wp,  0.5_wp,  0.5_wp,  0.0_wp,  0.0_wp, 0.0_wp, 0.0_wp  &
573                   /), (/maxspec,4/) )
574
575       REAL(wp), DIMENSION(1:3,1:4) ::  pm_frac_table = & !< rel. mass
576                                     RESHAPE( (/ 0.016_wp, 0.000_wp, 0.984_wp, &
577                                                 0.000_wp, 1.000_wp, 0.000_wp, &
578                                                 0.000_wp, 0.000_wp, 1.000_wp, &
579                                                 1.000_wp, 0.000_wp, 1.000_wp  &
580                                              /), (/3,4/) )
581
582    END TYPE salsa_emission_mode_type
583!
584!-- Aerosol emission data values
585    TYPE salsa_emission_value_type
586
587       REAL(wp) ::  fill  !< fill value
588
589       REAL(wp), DIMENSION(:), ALLOCATABLE :: preproc_mass_fracs  !< mass fractions
590
591       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: def_mass_fracs  !< mass fractions per emis. category
592
593       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: def_data      !< surface emission values in PM
594       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: preproc_data  !< surface emission values per bin
595
596    END TYPE salsa_emission_value_type
597!
598!-- Prognostic variable: Aerosol size bin information (number (#/m3) and mass (kg/m3) concentration)
599!-- and the concentration of gaseous tracers (#/m3). Gas tracers are contained sequentially in
600!-- dimension 4 as:
601!-- 1. H2SO4, 2. HNO3, 3. NH3, 4. OCNV (non-volatile organics), 5. OCSV (semi-volatile)
602    TYPE salsa_variable
603
604       REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  init  !<
605
606       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  diss_s     !<
607       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  flux_s     !<
608       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  source     !<
609       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  sums_ws_l  !<
610
611       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l  !<
612       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l  !<
613
614       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  conc     !<
615       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  conc_p   !<
616       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  tconc_m  !<
617
618    END TYPE salsa_variable
619!
620!-- Datatype used to store information about the binned size distributions of aerosols
621    TYPE t_section
622
623       REAL(wp) ::  dmid     !< bin middle diameter (m)
624       REAL(wp) ::  vhilim   !< bin volume at the high limit
625       REAL(wp) ::  vlolim   !< bin volume at the low limit
626       REAL(wp) ::  vratiohi !< volume ratio between the center and high limit
627       REAL(wp) ::  vratiolo !< volume ratio between the center and low limit
628       !******************************************************
629       ! ^ Do NOT change the stuff above after initialization !
630       !******************************************************
631       REAL(wp) ::  core    !< Volume of dry particle
632       REAL(wp) ::  dwet    !< Wet diameter or mean droplet diameter (m)
633       REAL(wp) ::  numc    !< Number concentration of particles/droplets (#/m3)
634       REAL(wp) ::  veqh2o  !< Equilibrium H2O concentration for each particle
635
636       REAL(wp), DIMENSION(maxspec+1) ::  volc !< Volume concentrations (m^3/m^3) of aerosols +
637                                               !< water. Since most of the stuff in SALSA is hard
638                                               !< coded, these *have to be* in the order
639                                               !< 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
640    END TYPE t_section
641
642    TYPE(salsa_emission_attribute_type) ::  aero_emission_att  !< emission attributes
643    TYPE(salsa_emission_value_type)     ::  aero_emission      !< emission values
644    TYPE(salsa_emission_mode_type)      ::  def_modes          !< default emission modes
645
646    TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
647    TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
648
649    TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
650
651    TYPE(match_surface) ::  lsm_to_depo_h  !< to match the deposition module and horizontal LSM surfaces
652    TYPE(match_surface) ::  usm_to_depo_h  !< to match the deposition module and horizontal USM surfaces
653
654    TYPE(match_surface), DIMENSION(0:3) ::  lsm_to_depo_v  !< to match the deposition mod. and vertical LSM surfaces
655    TYPE(match_surface), DIMENSION(0:3) ::  usm_to_depo_v  !< to match the deposition mod. and vertical USM surfaces
656!
657!-- SALSA variables: as x = x(k,j,i,bin).
658!-- The 4th dimension contains all the size bins sequentially for each aerosol species  + water.
659!
660!-- Prognostic variables:
661!
662!-- Number concentration (#/m3)
663    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_number  !<
664    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_1  !<
665    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_2  !<
666    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_3  !<
667!
668!-- Mass concentration (kg/m3)
669    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_mass  !<
670    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_1  !<
671    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_2  !<
672    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_3  !<
673!
674!-- Gaseous concentrations (#/m3)
675    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  salsa_gas  !<
676    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_1  !<
677    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_2  !<
678    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_3  !<
679!
680!-- Diagnostic tracers
681    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  sedim_vd  !< sedimentation velocity per bin (m/s)
682    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  ra_dry    !< aerosol dry radius (m)
683
684!-- Particle component index tables
685    TYPE(component_index) :: prtcl  !< Contains "getIndex" which gives the index for a given aerosol
686                                    !< component name: 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
687!
688!-- Data output arrays:
689!
690!-- Gases:
691    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_h2so4_av  !< H2SO4
692    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_hno3_av   !< HNO3
693    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_nh3_av    !< NH3
694    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_ocnv_av   !< non-volatile OC
695    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_ocsv_av   !< semi-volatile OC
696!
697!-- Integrated:
698    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ldsa_av  !< lung-deposited surface area
699    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ntot_av  !< total number concentration
700    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  pm25_av  !< PM2.5
701    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  pm10_av  !< PM10
702!
703!-- In the particle phase:
704    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_bc_av   !< black carbon
705    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_du_av   !< dust
706    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_h2o_av  !< liquid water
707    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_nh_av   !< ammonia
708    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_no_av   !< nitrates
709    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_oc_av   !< org. carbon
710    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_so4_av  !< sulphates
711    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_ss_av   !< sea salt
712!
713!-- Bin specific mass and number concentrations:
714    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mbins_av  !< bin mas
715    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nbins_av  !< bin number
716
717!
718!-- PALM interfaces:
719!
720!-- Boundary conditions:
721    INTERFACE salsa_boundary_conds
722       MODULE PROCEDURE salsa_boundary_conds
723       MODULE PROCEDURE salsa_boundary_conds_decycle
724    END INTERFACE salsa_boundary_conds
725!
726!-- Data output checks for 2D/3D data to be done in check_parameters
727    INTERFACE salsa_check_data_output
728       MODULE PROCEDURE salsa_check_data_output
729    END INTERFACE salsa_check_data_output
730!
731!-- Input parameter checks to be done in check_parameters
732    INTERFACE salsa_check_parameters
733       MODULE PROCEDURE salsa_check_parameters
734    END INTERFACE salsa_check_parameters
735!
736!-- Averaging of 3D data for output
737    INTERFACE salsa_3d_data_averaging
738       MODULE PROCEDURE salsa_3d_data_averaging
739    END INTERFACE salsa_3d_data_averaging
740!
741!-- Data output of 2D quantities
742    INTERFACE salsa_data_output_2d
743       MODULE PROCEDURE salsa_data_output_2d
744    END INTERFACE salsa_data_output_2d
745!
746!-- Data output of 3D data
747    INTERFACE salsa_data_output_3d
748       MODULE PROCEDURE salsa_data_output_3d
749    END INTERFACE salsa_data_output_3d
750!
751!-- Data output of 3D data
752    INTERFACE salsa_data_output_mask
753       MODULE PROCEDURE salsa_data_output_mask
754    END INTERFACE salsa_data_output_mask
755!
756!-- Definition of data output quantities
757    INTERFACE salsa_define_netcdf_grid
758       MODULE PROCEDURE salsa_define_netcdf_grid
759    END INTERFACE salsa_define_netcdf_grid
760!
761!-- Output of information to the header file
762    INTERFACE salsa_header
763       MODULE PROCEDURE salsa_header
764    END INTERFACE salsa_header
765!
766!-- Initialization actions
767    INTERFACE salsa_init
768       MODULE PROCEDURE salsa_init
769    END INTERFACE salsa_init
770!
771!-- Initialization of arrays
772    INTERFACE salsa_init_arrays
773       MODULE PROCEDURE salsa_init_arrays
774    END INTERFACE salsa_init_arrays
775!
776!-- Writing of binary output for restart runs  !!! renaming?!
777    INTERFACE salsa_wrd_local
778       MODULE PROCEDURE salsa_wrd_local
779    END INTERFACE salsa_wrd_local
780!
781!-- Reading of NAMELIST parameters
782    INTERFACE salsa_parin
783       MODULE PROCEDURE salsa_parin
784    END INTERFACE salsa_parin
785!
786!-- Reading of parameters for restart runs
787    INTERFACE salsa_rrd_local
788       MODULE PROCEDURE salsa_rrd_local
789    END INTERFACE salsa_rrd_local
790!
791!-- Swapping of time levels (required for prognostic variables)
792    INTERFACE salsa_swap_timelevel
793       MODULE PROCEDURE salsa_swap_timelevel
794    END INTERFACE salsa_swap_timelevel
795!
796!-- Interface between PALM and salsa
797    INTERFACE salsa_driver
798       MODULE PROCEDURE salsa_driver
799    END INTERFACE salsa_driver
800
801!-- Actions salsa variables
802    INTERFACE salsa_actions
803       MODULE PROCEDURE salsa_actions
804       MODULE PROCEDURE salsa_actions_ij
805    END INTERFACE salsa_actions
806!
807!-- Non-advective processes (i.e. aerosol dynamic reactions) for salsa variables
808    INTERFACE salsa_non_advective_processes
809       MODULE PROCEDURE salsa_non_advective_processes
810       MODULE PROCEDURE salsa_non_advective_processes_ij
811    END INTERFACE salsa_non_advective_processes
812!
813!-- Exchange horiz for salsa variables
814    INTERFACE salsa_exchange_horiz_bounds
815       MODULE PROCEDURE salsa_exchange_horiz_bounds
816    END INTERFACE salsa_exchange_horiz_bounds
817!
818!-- Prognostics equations for salsa variables
819    INTERFACE salsa_prognostic_equations
820       MODULE PROCEDURE salsa_prognostic_equations
821       MODULE PROCEDURE salsa_prognostic_equations_ij
822    END INTERFACE salsa_prognostic_equations
823!
824!-- Tendency salsa variables
825    INTERFACE salsa_tendency
826       MODULE PROCEDURE salsa_tendency
827       MODULE PROCEDURE salsa_tendency_ij
828    END INTERFACE salsa_tendency
829
830
831    SAVE
832
833    PRIVATE
834!
835!-- Public functions:
836    PUBLIC salsa_boundary_conds, salsa_check_data_output, salsa_check_parameters,                  &
837           salsa_3d_data_averaging, salsa_data_output_2d, salsa_data_output_3d,                    &
838           salsa_data_output_mask, salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver,      &
839           salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin,        &
840           salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local,     &
841           salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds
842!
843!-- Public parameters, constants and initial values
844    PUBLIC bc_am_t_val, bc_an_t_val, bc_gt_t_val, dots_salsa, dt_salsa,                            &
845           ibc_salsa_b, last_salsa_time, lsdepo, nest_salsa, salsa, salsa_gases_from_chem,         &
846           skip_time_do_salsa
847!
848!-- Public prognostic variables
849    PUBLIC aerosol_mass, aerosol_number, gconc_2, mconc_2, nbins_aerosol, ncc, ncomponents_mass,   &
850           nclim, nconc_2, ngases_salsa, prtcl, ra_dry, salsa_gas, sedim_vd
851
852
853 CONTAINS
854
855!------------------------------------------------------------------------------!
856! Description:
857! ------------
858!> Parin for &salsa_par for new modules
859!------------------------------------------------------------------------------!
860 SUBROUTINE salsa_parin
861
862    IMPLICIT NONE
863
864    CHARACTER(LEN=80) ::  line   !< dummy string that contains the current line
865                                  !< of the parameter file
866
867    NAMELIST /salsa_parameters/      aerosol_flux_dpg, aerosol_flux_mass_fracs_a,                  &
868                                     aerosol_flux_mass_fracs_b, aerosol_flux_sigmag,               &
869                                     advect_particle_water, bc_salsa_b, bc_salsa_t, decycle_lr,    &
870                                     decycle_method, decycle_ns, depo_pcm_par, depo_pcm_type,      &
871                                     depo_surf_par, dpg, dt_salsa, feedback_to_palm, h2so4_init,   &
872                                     hno3_init, init_gases_type, init_aerosol_type, listspec,      &
873                                     mass_fracs_a, mass_fracs_b, n_lognorm, nbin, nest_salsa, nf2a,&
874                                     nh3_init, nj3, nlcnd, nlcndgas, nlcndh2oae, nlcoag, nldepo,   &
875                                     nldepo_pcm,  nldepo_surf, nldistupdate, nsnucl, ocnv_init,    &
876                                     ocsv_init, read_restart_data_salsa, reglim, salsa,            &
877                                     salsa_emission_mode, sigmag, skip_time_do_salsa,              &
878                                     surface_aerosol_flux, van_der_waals_coagc, write_binary_salsa
879
880    line = ' '
881!
882!-- Try to find salsa package
883    REWIND ( 11 )
884    line = ' '
885    DO WHILE ( INDEX( line, '&salsa_parameters' ) == 0 )
886       READ ( 11, '(A)', END=10 )  line
887    ENDDO
888    BACKSPACE ( 11 )
889!
890!-- Read user-defined namelist
891    READ ( 11, salsa_parameters )
892!
893!-- Enable salsa (salsa switch in modules.f90)
894    salsa = .TRUE.
895
896 10 CONTINUE
897
898 END SUBROUTINE salsa_parin
899
900!------------------------------------------------------------------------------!
901! Description:
902! ------------
903!> Check parameters routine for salsa.
904!------------------------------------------------------------------------------!
905 SUBROUTINE salsa_check_parameters
906
907    USE control_parameters,                                                                        &
908        ONLY:  message_string
909
910    IMPLICIT NONE
911
912!
913!-- Checks go here (cf. check_parameters.f90).
914    IF ( salsa  .AND.  .NOT.  humidity )  THEN
915       WRITE( message_string, * ) 'salsa = ', salsa, ' is not allowed with humidity = ', humidity
916       CALL message( 'salsa_check_parameters', 'PA0594', 1, 2, 0, 6, 0 )
917    ENDIF
918
919    IF ( bc_salsa_b == 'dirichlet' )  THEN
920       ibc_salsa_b = 0
921    ELSEIF ( bc_salsa_b == 'neumann' )  THEN
922       ibc_salsa_b = 1
923    ELSE
924       message_string = 'unknown boundary condition: bc_salsa_b = "' // TRIM( bc_salsa_t ) // '"'
925       CALL message( 'salsa_check_parameters', 'PA0595', 1, 2, 0, 6, 0 )
926    ENDIF
927
928    IF ( bc_salsa_t == 'dirichlet' )  THEN
929       ibc_salsa_t = 0
930    ELSEIF ( bc_salsa_t == 'neumann' )  THEN
931       ibc_salsa_t = 1
932    ELSEIF ( bc_salsa_t == 'nested' )  THEN
933       ibc_salsa_t = 2
934    ELSE
935       message_string = 'unknown boundary condition: bc_salsa_t = "' // TRIM( bc_salsa_t ) // '"'
936       CALL message( 'salsa_check_parameters', 'PA0596', 1, 2, 0, 6, 0 )
937    ENDIF
938
939    IF ( nj3 < 1  .OR.  nj3 > 3 )  THEN
940       message_string = 'unknown nj3 (must be 1-3)'
941       CALL message( 'salsa_check_parameters', 'PA0597', 1, 2, 0, 6, 0 )
942    ENDIF
943
944    IF ( salsa_emission_mode /= 'no_emission'  .AND.  ibc_salsa_b  == 0 ) THEN
945       message_string = 'salsa_emission_mode /= "no_emission" requires bc_salsa_b = "Neumann"'
946       CALL message( 'salsa_check_parameters','PA0598', 1, 2, 0, 6, 0 )
947    ENDIF
948
949    IF ( salsa_emission_mode /= 'no_emission' )  include_emission = .TRUE.
950
951 END SUBROUTINE salsa_check_parameters
952
953!------------------------------------------------------------------------------!
954!
955! Description:
956! ------------
957!> Subroutine defining appropriate grid for netcdf variables.
958!> It is called out from subroutine netcdf.
959!> Same grid as for other scalars (see netcdf_interface_mod.f90)
960!------------------------------------------------------------------------------!
961 SUBROUTINE salsa_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
962
963    IMPLICIT NONE
964
965    CHARACTER(LEN=*), INTENT(OUT) ::  grid_x   !<
966    CHARACTER(LEN=*), INTENT(OUT) ::  grid_y   !<
967    CHARACTER(LEN=*), INTENT(OUT) ::  grid_z   !<
968    CHARACTER(LEN=*), INTENT(IN)  ::  var      !<
969
970    LOGICAL, INTENT(OUT) ::  found   !<
971
972    found  = .TRUE.
973!
974!-- Check for the grid
975
976    IF ( var(1:2) == 'g_' )  THEN
977       grid_x = 'x'
978       grid_y = 'y'
979       grid_z = 'zu'
980    ELSEIF ( var(1:4) == 'LDSA' )  THEN
981       grid_x = 'x'
982       grid_y = 'y'
983       grid_z = 'zu'
984    ELSEIF ( var(1:5) == 'm_bin' )  THEN
985       grid_x = 'x'
986       grid_y = 'y'
987       grid_z = 'zu'
988    ELSEIF ( var(1:5) == 'N_bin' )  THEN
989       grid_x = 'x'
990       grid_y = 'y'
991       grid_z = 'zu'
992    ELSEIF ( var(1:4) == 'Ntot' ) THEN
993       grid_x = 'x'
994       grid_y = 'y'
995       grid_z = 'zu'
996    ELSEIF ( var(1:2) == 'PM' )  THEN
997       grid_x = 'x'
998       grid_y = 'y'
999       grid_z = 'zu'
1000    ELSEIF ( var(1:2) == 's_' )  THEN
1001       grid_x = 'x'
1002       grid_y = 'y'
1003       grid_z = 'zu'
1004    ELSE
1005       found  = .FALSE.
1006       grid_x = 'none'
1007       grid_y = 'none'
1008       grid_z = 'none'
1009    ENDIF
1010
1011 END SUBROUTINE salsa_define_netcdf_grid
1012
1013!------------------------------------------------------------------------------!
1014! Description:
1015! ------------
1016!> Header output for new module
1017!------------------------------------------------------------------------------!
1018 SUBROUTINE salsa_header( io )
1019
1020    USE indices,                                                                                   &
1021        ONLY:  nx, ny, nz
1022
1023    IMPLICIT NONE
1024 
1025    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
1026!
1027!-- Write SALSA header
1028    WRITE( io, 1 )
1029    WRITE( io, 2 ) skip_time_do_salsa
1030    WRITE( io, 3 ) dt_salsa
1031    WRITE( io, 4 )  nz, ny, nx, nbins_aerosol
1032    IF ( advect_particle_water )  THEN
1033       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol,              &
1034                        advect_particle_water
1035    ELSE
1036       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncc*nbins_aerosol, advect_particle_water
1037    ENDIF
1038    IF ( .NOT. salsa_gases_from_chem )  THEN
1039       WRITE( io, 6 )  SHAPE( aerosol_mass(1)%conc ), ngases_salsa, salsa_gases_from_chem
1040    ENDIF
1041    WRITE( io, 7 )
1042    IF ( nsnucl > 0 )   WRITE( io, 8 ) nsnucl, nj3
1043    IF ( nlcoag )       WRITE( io, 9 )
1044    IF ( nlcnd )        WRITE( io, 10 ) nlcndgas, nlcndh2oae
1045    IF ( lspartition )  WRITE( io, 11 )
1046    IF ( nldepo )       WRITE( io, 12 ) nldepo_pcm, nldepo_surf
1047    WRITE( io, 13 )  reglim, nbin, bin_low_limits
1048    IF ( init_aerosol_type == 0 )  WRITE( io, 14 ) nsect
1049    WRITE( io, 15 ) ncc, listspec, mass_fracs_a, mass_fracs_b
1050    IF ( .NOT. salsa_gases_from_chem )  THEN
1051       WRITE( io, 16 ) ngases_salsa, h2so4_init, hno3_init, nh3_init, ocnv_init, ocsv_init
1052    ENDIF
1053    WRITE( io, 17 )  init_aerosol_type, init_gases_type
1054    IF ( init_aerosol_type == 0 )  THEN
1055       WRITE( io, 18 )  dpg, sigmag, n_lognorm
1056    ELSE
1057       WRITE( io, 19 )
1058    ENDIF
1059    IF ( nest_salsa )  WRITE( io, 20 )  nest_salsa
1060    WRITE( io, 21 ) salsa_emission_mode
1061    IF ( salsa_emission_mode == 'uniform' )  THEN
1062       WRITE( io, 22 ) surface_aerosol_flux, aerosol_flux_dpg, aerosol_flux_sigmag,                &
1063                       aerosol_flux_mass_fracs_a
1064    ENDIF
1065    IF ( SUM( aerosol_flux_mass_fracs_b ) > 0.0_wp  .OR. salsa_emission_mode == 'read_from_file' ) &
1066    THEN
1067       WRITE( io, 23 )
1068    ENDIF
1069
10701   FORMAT (//' SALSA information:'/                                                               &
1071              ' ------------------------------'/)
10722   FORMAT   ('    Starts at: skip_time_do_salsa = ', F10.2, '  s')
10733   FORMAT  (/'    Timestep: dt_salsa = ', F6.2, '  s')
10744   FORMAT  (/'    Array shape (z,y,x,bins):'/                                                     &
1075              '       aerosol_number:  ', 4(I3)) 
10765   FORMAT  (/'       aerosol_mass:    ', 4(I3),/                                                  &
1077              '       (advect_particle_water = ', L1, ')')
10786   FORMAT   ('       salsa_gas: ', 4(I3),/                                                        &
1079              '       (salsa_gases_from_chem = ', L1, ')')
10807   FORMAT  (/'    Aerosol dynamic processes included: ')
10818   FORMAT  (/'       nucleation (scheme = ', I1, ' and J3 parametrization = ', I1, ')')
10829   FORMAT  (/'       coagulation')
108310  FORMAT  (/'       condensation (of precursor gases = ', L1, ' and water vapour = ', L1, ')' )
108411  FORMAT  (/'       dissolutional growth by HNO3 and NH3')
108512  FORMAT  (/'       dry deposition (on vegetation = ', L1, ' and on topography = ', L1, ')')
108613  FORMAT  (/'    Aerosol bin subrange limits (in metres): ',  3(ES10.2E3), /                     &
1087              '    Number of size bins for each aerosol subrange: ', 2I3,/                         &
1088              '    Aerosol bin limits (in metres): ', 9(ES10.2E3))
108914  FORMAT   ('    Initial number concentration in bins at the lowest level (#/m**3):', 9(ES10.2E3))
109015  FORMAT  (/'    Number of chemical components used: ', I1,/                                     &
1091              '       Species: ',7(A6),/                                                           &
1092              '    Initial relative contribution of each species to particle volume in:',/         &
1093              '       a-bins: ', 7(F6.3),/                                                         &
1094              '       b-bins: ', 7(F6.3))
109516  FORMAT  (/'    Number of gaseous tracers used: ', I1,/                                         &
1096              '    Initial gas concentrations:',/                                                  &
1097              '       H2SO4: ',ES12.4E3, ' #/m**3',/                                               &
1098              '       HNO3:  ',ES12.4E3, ' #/m**3',/                                               &
1099              '       NH3:   ',ES12.4E3, ' #/m**3',/                                               &
1100              '       OCNV:  ',ES12.4E3, ' #/m**3',/                                               &
1101              '       OCSV:  ',ES12.4E3, ' #/m**3')
110217   FORMAT (/'   Initialising concentrations: ', /                                                &
1103              '      Aerosol size distribution: init_aerosol_type = ', I1,/                        &
1104              '      Gas concentrations: init_gases_type = ', I1 )
110518   FORMAT ( '      Mode diametres: dpg(nmod) = ', 7(F7.3), ' (m)', /                             &
1106              '      Standard deviation: sigmag(nmod) = ', 7(F7.2),/                               &
1107              '      Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3), ' (#/m3)' )
110819   FORMAT (/'      Size distribution read from a file.')
110920   FORMAT (/'   Nesting for salsa variables: ', L1 )
111021   FORMAT (/'   Emissions: salsa_emission_mode = ', A )
111122   FORMAT (/'      surface_aerosol_flux = ', ES12.4E3, ' #/m**2/s', /                            &
1112              '      aerosol_flux_dpg     =  ', 7(F7.3), ' (m)', /                                 &
1113              '      aerosol_flux_sigmag  =  ', 7(F7.2), /                                         &
1114              '      aerosol_mass_fracs_a =  ', 7(ES12.4E3) )
111523   FORMAT (/'      (currently all emissions are soluble!)')
1116
1117 END SUBROUTINE salsa_header
1118
1119!------------------------------------------------------------------------------!
1120! Description:
1121! ------------
1122!> Allocate SALSA arrays and define pointers if required
1123!------------------------------------------------------------------------------!
1124 SUBROUTINE salsa_init_arrays
1125
1126    USE chem_gasphase_mod,                                                                         &
1127        ONLY:  nvar
1128
1129    USE surface_mod,                                                                               &
1130        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
1131
1132    IMPLICIT NONE
1133
1134    INTEGER(iwp) ::  gases_available !< Number of available gas components in the chemistry model
1135    INTEGER(iwp) ::  i               !< loop index for allocating
1136    INTEGER(iwp) ::  l               !< loop index for allocating: surfaces
1137    INTEGER(iwp) ::  lsp             !< loop index for chem species in the chemistry model
1138
1139    gases_available = 0
1140!
1141!-- Allocate prognostic variables (see salsa_swap_timelevel)
1142!
1143!-- Set derived indices:
1144!-- (This does the same as the subroutine salsa_initialize in SALSA/UCLALES-SALSA)
1145    start_subrange_1a = 1  ! 1st index of subrange 1a
1146    start_subrange_2a = start_subrange_1a + nbin(1)  ! 1st index of subrange 2a
1147    end_subrange_1a   = start_subrange_2a - 1        ! last index of subrange 1a
1148    end_subrange_2a   = end_subrange_1a + nbin(2)    ! last index of subrange 2a
1149
1150!
1151!-- If the fraction of insoluble aerosols in subrange 2 is zero: do not allocate arrays for them
1152    IF ( nf2a > 0.999999_wp  .AND.  SUM( mass_fracs_b ) < 0.00001_wp )  THEN
1153       no_insoluble = .TRUE.
1154       start_subrange_2b = end_subrange_2a+1  ! 1st index of subrange 2b
1155       end_subrange_2b   = end_subrange_2a    ! last index of subrange 2b
1156    ELSE
1157       start_subrange_2b = start_subrange_2a + nbin(2)  ! 1st index of subrange 2b
1158       end_subrange_2b   = end_subrange_2a + nbin(2)    ! last index of subrange 2b
1159    ENDIF
1160
1161    nbins_aerosol = end_subrange_2b   ! total number of aerosol size bins
1162!
1163!-- Create index tables for different aerosol components
1164    CALL component_index_constructor( prtcl, ncc, maxspec, listspec )
1165
1166    ncomponents_mass = ncc
1167    IF ( advect_particle_water )  ncomponents_mass = ncc + 1  ! Add water
1168
1169!
1170!-- Allocate:
1171    ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass),                    &
1172              bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),&
1173              nsect(nbins_aerosol), massacc(nbins_aerosol) )
1174    ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) )
1175    IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1176    ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1177
1178!
1179!-- Aerosol number concentration
1180    ALLOCATE( aerosol_number(nbins_aerosol) )
1181    ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
1182              nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
1183              nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1184    nconc_1 = 0.0_wp
1185    nconc_2 = 0.0_wp
1186    nconc_3 = 0.0_wp
1187
1188    DO i = 1, nbins_aerosol
1189       aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => nconc_1(:,:,:,i)
1190       aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => nconc_2(:,:,:,i)
1191       aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i)
1192       ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                         &
1193                 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                         &
1194                 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
1195                 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
1196                 aerosol_number(i)%init(nzb:nzt+1),                                                &
1197                 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1198       aerosol_number(i)%init = nclim
1199       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
1200          ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
1201       ENDIF
1202    ENDDO
1203
1204!
1205!-- Aerosol mass concentration
1206    ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) )
1207    ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
1208              mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
1209              mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) )
1210    mconc_1 = 0.0_wp
1211    mconc_2 = 0.0_wp
1212    mconc_3 = 0.0_wp
1213
1214    DO i = 1, ncomponents_mass*nbins_aerosol
1215       aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => mconc_1(:,:,:,i)
1216       aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => mconc_2(:,:,:,i)
1217       aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i)
1218       ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                           &
1219                 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                           &
1220                 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
1221                 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
1222                 aerosol_mass(i)%init(nzb:nzt+1),                                                  &
1223                 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
1224       aerosol_mass(i)%init = mclim
1225       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
1226          ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
1227       ENDIF
1228    ENDDO
1229
1230!
1231!-- Surface fluxes: answs = aerosol number, amsws = aerosol mass
1232!
1233!-- Horizontal surfaces: default type
1234    DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
1235       ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) )
1236       ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) )
1237       surf_def_h(l)%answs = 0.0_wp
1238       surf_def_h(l)%amsws = 0.0_wp
1239    ENDDO
1240!
1241!-- Horizontal surfaces: natural type
1242    ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) )
1243    ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) )
1244    surf_lsm_h%answs = 0.0_wp
1245    surf_lsm_h%amsws = 0.0_wp
1246!
1247!-- Horizontal surfaces: urban type
1248    ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) )
1249    ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) )
1250    surf_usm_h%answs = 0.0_wp
1251    surf_usm_h%amsws = 0.0_wp
1252
1253!
1254!-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing
1255    DO  l = 0, 3
1256       ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) )
1257       surf_def_v(l)%answs = 0.0_wp
1258       ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1259       surf_def_v(l)%amsws = 0.0_wp
1260
1261       ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) )
1262       surf_lsm_v(l)%answs = 0.0_wp
1263       ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1264       surf_lsm_v(l)%amsws = 0.0_wp
1265
1266       ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) )
1267       surf_usm_v(l)%answs = 0.0_wp
1268       ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1269       surf_usm_v(l)%amsws = 0.0_wp
1270
1271    ENDDO
1272
1273!
1274!-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV)
1275!-- (number concentration (#/m3) )
1276!
1277!-- If chemistry is on, read gas phase concentrations from there. Otherwise,
1278!-- allocate salsa_gas array.
1279
1280    IF ( air_chemistry )  THEN
1281       DO  lsp = 1, nvar
1282          SELECT CASE ( TRIM( chem_species(lsp)%name ) )
1283             CASE ( 'H2SO4', 'h2so4' )
1284                gases_available = gases_available + 1
1285                gas_index_chem(1) = lsp
1286             CASE ( 'HNO3', 'hno3' )
1287                gases_available = gases_available + 1
1288                gas_index_chem(2) = lsp
1289             CASE ( 'NH3', 'nh3' )
1290                gases_available = gases_available + 1
1291                gas_index_chem(3) = lsp
1292             CASE ( 'OCNV', 'ocnv' )
1293                gases_available = gases_available + 1
1294                gas_index_chem(4) = lsp
1295             CASE ( 'OCSV', 'ocsv' )
1296                gases_available = gases_available + 1
1297                gas_index_chem(5) = lsp
1298          END SELECT
1299       ENDDO
1300
1301       IF ( gases_available == ngases_salsa )  THEN
1302          salsa_gases_from_chem = .TRUE.
1303       ELSE
1304          WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// &
1305                                     'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)'
1306       CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 )
1307       ENDIF
1308
1309    ELSE
1310
1311       ALLOCATE( salsa_gas(ngases_salsa) )
1312       ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
1313                 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
1314                 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
1315       gconc_1 = 0.0_wp
1316       gconc_2 = 0.0_wp
1317       gconc_3 = 0.0_wp
1318
1319       DO i = 1, ngases_salsa
1320          salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => gconc_1(:,:,:,i)
1321          salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => gconc_2(:,:,:,i)
1322          salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i)
1323          ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
1324                    salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
1325                    salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1326                    salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1327                    salsa_gas(i)%init(nzb:nzt+1),                              &
1328                    salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1329          salsa_gas(i)%init = nclim
1330          IF ( include_emission )  ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
1331       ENDDO
1332!
1333!--    Surface fluxes: gtsws = gaseous tracer flux
1334!
1335!--    Horizontal surfaces: default type
1336       DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
1337          ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) )
1338          surf_def_h(l)%gtsws = 0.0_wp
1339       ENDDO
1340!--    Horizontal surfaces: natural type
1341       ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) )
1342       surf_lsm_h%gtsws = 0.0_wp
1343!--    Horizontal surfaces: urban type
1344       ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) )
1345       surf_usm_h%gtsws = 0.0_wp
1346!
1347!--    Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
1348!--    westward (l=3) facing
1349       DO  l = 0, 3
1350          ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) )
1351          surf_def_v(l)%gtsws = 0.0_wp
1352          ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) )
1353          surf_lsm_v(l)%gtsws = 0.0_wp
1354          ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) )
1355          surf_usm_v(l)%gtsws = 0.0_wp
1356       ENDDO
1357    ENDIF
1358
1359    IF ( ws_scheme_sca )  THEN
1360
1361       IF ( salsa )  THEN
1362          ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1363          sums_salsa_ws_l = 0.0_wp
1364       ENDIF
1365
1366    ENDIF
1367
1368 END SUBROUTINE salsa_init_arrays
1369
1370!------------------------------------------------------------------------------!
1371! Description:
1372! ------------
1373!> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA.
1374!> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are
1375!> also merged here.
1376!------------------------------------------------------------------------------!
1377 SUBROUTINE salsa_init
1378
1379    IMPLICIT NONE
1380
1381    INTEGER(iwp) :: i   !<
1382    INTEGER(iwp) :: ib  !< loop index for aerosol number bins
1383    INTEGER(iwp) :: ic  !< loop index for aerosol mass bins
1384    INTEGER(iwp) :: ig  !< loop index for gases
1385    INTEGER(iwp) :: ii  !< index for indexing
1386    INTEGER(iwp) :: j   !<
1387
1388    IF ( debug_output )  CALL debug_message( 'salsa_init', 'start' )
1389
1390    bin_low_limits = 0.0_wp
1391    k_topo_top     = 0
1392    nsect          = 0.0_wp
1393    massacc        = 1.0_wp
1394
1395!
1396!-- Indices for chemical components used (-1 = not used)
1397    ii = 0
1398    IF ( is_used( prtcl, 'SO4' ) )  THEN
1399       index_so4 = get_index( prtcl,'SO4' )
1400       ii = ii + 1
1401    ENDIF
1402    IF ( is_used( prtcl,'OC' ) )  THEN
1403       index_oc = get_index(prtcl, 'OC')
1404       ii = ii + 1
1405    ENDIF
1406    IF ( is_used( prtcl, 'BC' ) )  THEN
1407       index_bc = get_index( prtcl, 'BC' )
1408       ii = ii + 1
1409    ENDIF
1410    IF ( is_used( prtcl, 'DU' ) )  THEN
1411       index_du = get_index( prtcl, 'DU' )
1412       ii = ii + 1
1413    ENDIF
1414    IF ( is_used( prtcl, 'SS' ) )  THEN
1415       index_ss = get_index( prtcl, 'SS' )
1416       ii = ii + 1
1417    ENDIF
1418    IF ( is_used( prtcl, 'NO' ) )  THEN
1419       index_no = get_index( prtcl, 'NO' )
1420       ii = ii + 1
1421    ENDIF
1422    IF ( is_used( prtcl, 'NH' ) )  THEN
1423       index_nh = get_index( prtcl, 'NH' )
1424       ii = ii + 1
1425    ENDIF
1426!
1427!-- All species must be known
1428    IF ( ii /= ncc )  THEN
1429       message_string = 'Unknown aerosol species/component(s) given in the initialization'
1430       CALL message( 'salsa_mod: salsa_init', 'PA0600', 1, 2, 0, 6, 0 )
1431    ENDIF
1432!
1433!-- Partition and dissolutional growth by gaseous HNO3 and NH3
1434    IF ( index_no > 0  .AND.  index_nh > 0  .AND.  index_so4 > 0 )  lspartition = .TRUE.
1435!
1436!-- Initialise
1437!
1438!-- Aerosol size distribution (TYPE t_section)
1439    aero(:)%dwet     = 1.0E-10_wp
1440    aero(:)%veqh2o   = 1.0E-10_wp
1441    aero(:)%numc     = nclim
1442    aero(:)%core     = 1.0E-10_wp
1443    DO ic = 1, maxspec+1    ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
1444       aero(:)%volc(ic) = 0.0_wp
1445    ENDDO
1446
1447    IF ( nldepo )  sedim_vd = 0.0_wp
1448
1449    IF ( .NOT. salsa_gases_from_chem )  THEN
1450       IF ( .NOT. read_restart_data_salsa )  THEN
1451          salsa_gas(1)%conc = h2so4_init
1452          salsa_gas(2)%conc = hno3_init
1453          salsa_gas(3)%conc = nh3_init
1454          salsa_gas(4)%conc = ocnv_init
1455          salsa_gas(5)%conc = ocsv_init
1456       ENDIF
1457       DO  ig = 1, ngases_salsa
1458          salsa_gas(ig)%conc_p    = 0.0_wp
1459          salsa_gas(ig)%tconc_m   = 0.0_wp
1460          salsa_gas(ig)%flux_s    = 0.0_wp
1461          salsa_gas(ig)%diss_s    = 0.0_wp
1462          salsa_gas(ig)%flux_l    = 0.0_wp
1463          salsa_gas(ig)%diss_l    = 0.0_wp
1464          salsa_gas(ig)%sums_ws_l = 0.0_wp
1465          salsa_gas(ig)%conc_p    = salsa_gas(ig)%conc
1466       ENDDO
1467!
1468!--    Set initial value for gas compound tracer
1469       salsa_gas(1)%init = h2so4_init
1470       salsa_gas(2)%init = hno3_init
1471       salsa_gas(3)%init = nh3_init
1472       salsa_gas(4)%init = ocnv_init
1473       salsa_gas(5)%init = ocsv_init
1474    ENDIF
1475!
1476!-- Aerosol radius in each bin: dry and wet (m)
1477    ra_dry = 1.0E-10_wp
1478!
1479!-- Initialise aerosol tracers
1480    aero(:)%vhilim   = 0.0_wp
1481    aero(:)%vlolim   = 0.0_wp
1482    aero(:)%vratiohi = 0.0_wp
1483    aero(:)%vratiolo = 0.0_wp
1484    aero(:)%dmid     = 0.0_wp
1485!
1486!-- Initialise the sectional particle size distribution
1487    CALL set_sizebins
1488!
1489!-- Initialise location-dependent aerosol size distributions and chemical compositions:
1490    CALL aerosol_init
1491
1492!-- Initalisation run of SALSA + calculate the vertical top index of the topography
1493    DO  i = nxl, nxr
1494       DO  j = nys, nyn
1495
1496          k_topo_top(j,i) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,j,i), 12 ) ), DIM = 1 ) - 1
1497
1498          CALL salsa_driver( i, j, 1 )
1499          CALL salsa_diagnostics( i, j )
1500       ENDDO
1501    ENDDO
1502
1503    DO  ib = 1, nbins_aerosol
1504       aerosol_number(ib)%conc_p    = aerosol_number(ib)%conc
1505       aerosol_number(ib)%tconc_m   = 0.0_wp
1506       aerosol_number(ib)%flux_s    = 0.0_wp
1507       aerosol_number(ib)%diss_s    = 0.0_wp
1508       aerosol_number(ib)%flux_l    = 0.0_wp
1509       aerosol_number(ib)%diss_l    = 0.0_wp
1510       aerosol_number(ib)%sums_ws_l = 0.0_wp
1511    ENDDO
1512    DO  ic = 1, ncomponents_mass*nbins_aerosol
1513       aerosol_mass(ic)%conc_p    = aerosol_mass(ic)%conc
1514       aerosol_mass(ic)%tconc_m   = 0.0_wp
1515       aerosol_mass(ic)%flux_s    = 0.0_wp
1516       aerosol_mass(ic)%diss_s    = 0.0_wp
1517       aerosol_mass(ic)%flux_l    = 0.0_wp
1518       aerosol_mass(ic)%diss_l    = 0.0_wp
1519       aerosol_mass(ic)%sums_ws_l = 0.0_wp
1520    ENDDO
1521!
1522!
1523!-- Initialise the deposition scheme and surface types
1524    IF ( nldepo )  CALL init_deposition
1525
1526    IF ( include_emission )  THEN
1527!
1528!--    Read in and initialize emissions
1529       CALL salsa_emission_setup( .TRUE. )
1530       IF ( .NOT. salsa_gases_from_chem  .AND.  salsa_emission_mode == 'read_from_file' )  THEN
1531          CALL salsa_gas_emission_setup( .TRUE. )
1532       ENDIF
1533    ENDIF
1534
1535    IF ( debug_output )  CALL debug_message( 'salsa_init', 'end' )
1536
1537 END SUBROUTINE salsa_init
1538
1539!------------------------------------------------------------------------------!
1540! Description:
1541! ------------
1542!> Initializes particle size distribution grid by calculating size bin limits
1543!> and mid-size for *dry* particles in each bin. Called from salsa_initialize
1544!> (only at the beginning of simulation).
1545!> Size distribution described using:
1546!>   1) moving center method (subranges 1 and 2)
1547!>      (Jacobson, Atmos. Env., 31, 131-144, 1997)
1548!>   2) fixed sectional method (subrange 3)
1549!> Size bins in each subrange are spaced logarithmically
1550!> based on given subrange size limits and bin number.
1551!
1552!> Mona changed 06/2017: Use geometric mean diameter to describe the mean
1553!> particle diameter in a size bin, not the arithmeric mean which clearly
1554!> overestimates the total particle volume concentration.
1555!
1556!> Coded by:
1557!> Hannele Korhonen (FMI) 2005
1558!> Harri Kokkola (FMI) 2006
1559!
1560!> Bug fixes for box model + updated for the new aerosol datatype:
1561!> Juha Tonttila (FMI) 2014
1562!------------------------------------------------------------------------------!
1563 SUBROUTINE set_sizebins
1564
1565    IMPLICIT NONE
1566
1567    INTEGER(iwp) ::  cc  !< running index
1568    INTEGER(iwp) ::  dd  !< running index
1569
1570    REAL(wp) ::  ratio_d  !< ratio of the upper and lower diameter of subranges
1571!
1572!-- vlolim&vhilim: min & max *dry* volumes [fxm]
1573!-- dmid: bin mid *dry* diameter (m)
1574!-- vratiolo&vratiohi: volume ratio between the center and low/high limit
1575!
1576!-- 1) Size subrange 1:
1577    ratio_d = reglim(2) / reglim(1)   ! section spacing (m)
1578    DO  cc = start_subrange_1a, end_subrange_1a
1579       aero(cc)%vlolim = api6 * ( reglim(1) * ratio_d**( REAL( cc-1 ) / nbin(1) ) )**3
1580       aero(cc)%vhilim = api6 * ( reglim(1) * ratio_d**( REAL( cc ) / nbin(1) ) )**3
1581       aero(cc)%dmid = SQRT( ( aero(cc)%vhilim / api6 )**0.33333333_wp *                           &
1582                             ( aero(cc)%vlolim / api6 )**0.33333333_wp )
1583       aero(cc)%vratiohi = aero(cc)%vhilim / ( api6 * aero(cc)%dmid**3 )
1584       aero(cc)%vratiolo = aero(cc)%vlolim / ( api6 * aero(cc)%dmid**3 )
1585    ENDDO
1586!
1587!-- 2) Size subrange 2:
1588!-- 2.1) Sub-subrange 2a: high hygroscopicity
1589    ratio_d = reglim(3) / reglim(2)   ! section spacing
1590    DO  dd = start_subrange_2a, end_subrange_2a
1591       cc = dd - start_subrange_2a
1592       aero(dd)%vlolim = api6 * ( reglim(2) * ratio_d**( REAL( cc ) / nbin(2) ) )**3
1593       aero(dd)%vhilim = api6 * ( reglim(2) * ratio_d**( REAL( cc+1 ) / nbin(2) ) )**3
1594       aero(dd)%dmid = SQRT( ( aero(dd)%vhilim / api6 )**0.33333333_wp *                           &
1595                             ( aero(dd)%vlolim / api6 )**0.33333333_wp )
1596       aero(dd)%vratiohi = aero(dd)%vhilim / ( api6 * aero(dd)%dmid**3 )
1597       aero(dd)%vratiolo = aero(dd)%vlolim / ( api6 * aero(dd)%dmid**3 )
1598    ENDDO
1599!
1600!-- 2.2) Sub-subrange 2b: low hygroscopicity
1601    IF ( .NOT. no_insoluble )  THEN
1602       aero(start_subrange_2b:end_subrange_2b)%vlolim   = aero(start_subrange_2a:end_subrange_2a)%vlolim
1603       aero(start_subrange_2b:end_subrange_2b)%vhilim   = aero(start_subrange_2a:end_subrange_2a)%vhilim
1604       aero(start_subrange_2b:end_subrange_2b)%dmid     = aero(start_subrange_2a:end_subrange_2a)%dmid
1605       aero(start_subrange_2b:end_subrange_2b)%vratiohi = aero(start_subrange_2a:end_subrange_2a)%vratiohi
1606       aero(start_subrange_2b:end_subrange_2b)%vratiolo = aero(start_subrange_2a:end_subrange_2a)%vratiolo
1607    ENDIF
1608!
1609!-- Initialize the wet diameter with the bin dry diameter to avoid numerical problems later
1610    aero(:)%dwet = aero(:)%dmid
1611!
1612!-- Save bin limits (lower diameter) to be delivered to PALM if needed
1613    DO cc = 1, nbins_aerosol
1614       bin_low_limits(cc) = ( aero(cc)%vlolim / api6 )**0.33333333_wp
1615    ENDDO
1616
1617 END SUBROUTINE set_sizebins
1618
1619!------------------------------------------------------------------------------!
1620! Description:
1621! ------------
1622!> Initilize altitude-dependent aerosol size distributions and compositions.
1623!>
1624!> Mona added 06/2017: Correct the number and mass concentrations by normalizing
1625!< by the given total number and mass concentration.
1626!>
1627!> Tomi Raatikainen, FMI, 29.2.2016
1628!------------------------------------------------------------------------------!
1629 SUBROUTINE aerosol_init
1630
1631    USE netcdf_data_input_mod,                                                                     &
1632        ONLY:  close_input_file, get_attribute, get_variable,                                      &
1633               netcdf_data_input_get_dimension_length, open_read_file
1634
1635    IMPLICIT NONE
1636
1637    CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name  !< chemical component name
1638
1639    INTEGER(iwp) ::  ee        !< index: end
1640    INTEGER(iwp) ::  i         !< loop index: x-direction
1641    INTEGER(iwp) ::  ib        !< loop index: size bins
1642    INTEGER(iwp) ::  ic        !< loop index: chemical components
1643    INTEGER(iwp) ::  id_dyn    !< NetCDF id of PIDS_DYNAMIC_SALSA
1644    INTEGER(iwp) ::  ig        !< loop index: gases
1645    INTEGER(iwp) ::  j         !< loop index: y-direction
1646    INTEGER(iwp) ::  k         !< loop index: z-direction
1647    INTEGER(iwp) ::  lod_aero  !< level of detail of inital aerosol concentrations
1648    INTEGER(iwp) ::  pr_nbins  !< Number of aerosol size bins in file
1649    INTEGER(iwp) ::  pr_ncc    !< Number of aerosol chemical components in file
1650    INTEGER(iwp) ::  pr_nz     !< Number of vertical grid-points in file
1651    INTEGER(iwp) ::  prunmode  !< running mode of SALSA
1652    INTEGER(iwp) ::  ss        !< index: start
1653
1654    INTEGER(iwp), DIMENSION(maxspec) ::  cc_input_to_model
1655
1656    LOGICAL  ::  netcdf_extend = .FALSE. !< Flag: netcdf file exists
1657
1658    REAL(wp) ::  flag  !< flag to mask topography grid points
1659
1660    REAL(wp), DIMENSION(nbins_aerosol) ::  core   !< size of the bin mid aerosol particle
1661
1662    REAL(wp), DIMENSION(0:nz+1) ::  pnf2a   !< number fraction in 2a
1663    REAL(wp), DIMENSION(0:nz+1) ::  pmfoc1a !< mass fraction of OC in 1a
1664
1665    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol)   ::  pndist  !< size dist as a function of height (#/m3)
1666    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2a   !< mass distributions in subrange 2a
1667    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2b   !< mass distributions in subrange 2b
1668
1669    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_dmid  !< vertical profile of aerosol bin diameters
1670    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_z     !< z levels of profiles
1671
1672    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_a  !< mass fraction: a
1673    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_b  !< and b
1674
1675    cc_input_to_model = 0
1676    prunmode = 1
1677!
1678!-- Bin mean aerosol particle volume (m3)
1679    core(:) = 0.0_wp
1680    core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3
1681!
1682!-- Set concentrations to zero
1683    pndist(:,:)  = 0.0_wp
1684    pnf2a(:)     = nf2a
1685    pmf2a(:,:)   = 0.0_wp
1686    pmf2b(:,:)   = 0.0_wp
1687    pmfoc1a(:)   = 0.0_wp
1688
1689    IF ( init_aerosol_type == 1 )  THEN
1690!
1691!--    Read input profiles from PIDS_DYNAMIC_SALSA
1692#if defined( __netcdf )
1693!
1694!--    Location-dependent size distributions and compositions.
1695       INQUIRE( FILE = TRIM( input_file_dynamic ) //  TRIM( coupling_char ), EXIST = netcdf_extend )
1696       IF ( netcdf_extend )  THEN
1697!
1698!--       Open file in read-only mode
1699          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
1700!
1701!--       Inquire dimensions:
1702          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' )
1703          IF ( pr_nz /= nz )  THEN
1704             WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//&
1705                                        'the number of numeric grid points.'
1706             CALL message( 'aerosol_init', 'PA0601', 1, 2, 0, 6, 0 )
1707          ENDIF
1708          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_ncc, 'composition_index' )
1709!
1710!--       Allocate memory
1711          ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                              &
1712                    pr_mass_fracs_b(nzb:nzt+1,pr_ncc) )
1713          pr_mass_fracs_a = 0.0_wp
1714          pr_mass_fracs_b = 0.0_wp
1715!
1716!--       Read vertical levels
1717          CALL get_variable( id_dyn, 'z', pr_z )
1718!
1719!--       Read name and index of chemical components
1720          CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc )
1721          DO  ic = 1, pr_ncc
1722             SELECT CASE ( TRIM( cc_name(ic) ) )
1723                CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' )
1724                   cc_input_to_model(1) = ic
1725                CASE ( 'OC', 'oc' )
1726                   cc_input_to_model(2) = ic
1727                CASE ( 'BC', 'bc' )
1728                   cc_input_to_model(3) = ic
1729                CASE ( 'DU', 'du' )
1730                   cc_input_to_model(4) = ic
1731                CASE ( 'SS', 'ss' )
1732                   cc_input_to_model(5) = ic
1733                CASE ( 'HNO3', 'hno3', 'NO', 'no' )
1734                   cc_input_to_model(6) = ic
1735                CASE ( 'NH3', 'nh3', 'NH', 'nh' )
1736                   cc_input_to_model(7) = ic
1737             END SELECT
1738          ENDDO
1739
1740          IF ( SUM( cc_input_to_model ) == 0 )  THEN
1741             message_string = 'None of the aerosol chemical components in ' // TRIM(               &
1742                              input_file_dynamic ) // ' correspond to ones applied in SALSA.'
1743             CALL message( 'salsa_mod: aerosol_init', 'PA0602', 2, 2, 0, 6, 0 )
1744          ENDIF
1745!
1746!--       Vertical profiles of mass fractions of different chemical components:
1747          CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a,              &
1748                             0, pr_ncc-1, 0, pr_nz-1 )
1749          CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_b', pr_mass_fracs_b,              &
1750                             0, pr_ncc-1, 0, pr_nz-1  )
1751!
1752!--       Match the input data with the chemical composition applied in the model
1753          DO  ic = 1, maxspec
1754             ss = cc_input_to_model(ic)
1755             IF ( ss == 0 )  CYCLE
1756             pmf2a(nzb+1:nzt+1,ic) = pr_mass_fracs_a(nzb:nzt,ss)
1757             pmf2b(nzb+1:nzt+1,ic) = pr_mass_fracs_b(nzb:nzt,ss)
1758          ENDDO
1759!
1760!--       Aerosol concentrations: lod=1 (total PM) or lod=2 (sectional number size distribution)
1761          CALL get_attribute( id_dyn, 'lod', lod_aero, .FALSE., 'init_atmosphere_aerosol' )
1762          IF ( lod_aero /= 1 )  THEN
1763             message_string = 'Currently only lod=1 accepted for init_atmosphere_aerosol'
1764             CALL message( 'salsa_mod: aerosol_init', 'PA0603', 2, 2, 0, 6, 0 )
1765          ELSE
1766!
1767!--          Bin mean diameters in the input file
1768             CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nbins, 'Dmid')
1769             IF ( pr_nbins /= nbins_aerosol )  THEN
1770                message_string = 'Number of size bins in init_atmosphere_aerosol does not match '  &
1771                                 // 'with that applied in the model'
1772                CALL message( 'salsa_mod: aerosol_init', 'PA0604', 2, 2, 0, 6, 0 )
1773             ENDIF
1774
1775             ALLOCATE( pr_dmid(pr_nbins) )
1776             pr_dmid    = 0.0_wp
1777
1778             CALL get_variable( id_dyn, 'Dmid', pr_dmid )
1779!
1780!--          Check whether the sectional representation conform to the one
1781!--          applied in the model
1782             IF ( ANY( ABS( ( aero(1:nbins_aerosol)%dmid - pr_dmid ) /                             &
1783                              aero(1:nbins_aerosol)%dmid )  > 0.1_wp )  ) THEN
1784                message_string = 'Mean diameters of the aerosol size bins in ' // TRIM(            &
1785                                 input_file_dynamic ) // ' do not match with the sectional '//     &
1786                                 'representation of the model.'
1787                CALL message( 'salsa_mod: aerosol_init', 'PA0605', 2, 2, 0, 6, 0 )
1788             ENDIF
1789!
1790!--          Inital aerosol concentrations
1791             CALL get_variable( id_dyn, 'init_atmosphere_aerosol', pndist(nzb+1:nzt,:),            &
1792                                0, pr_nbins-1, 0, pr_nz-1 )
1793          ENDIF
1794!
1795!--       Set bottom and top boundary condition (Neumann)
1796          pmf2a(nzb,:)    = pmf2a(nzb+1,:)
1797          pmf2a(nzt+1,:)  = pmf2a(nzt,:)
1798          pmf2b(nzb,:)    = pmf2b(nzb+1,:)
1799          pmf2b(nzt+1,:)  = pmf2b(nzt,:)
1800          pndist(nzb,:)   = pndist(nzb+1,:)
1801          pndist(nzt+1,:) = pndist(nzt,:)
1802
1803          IF ( index_so4 < 0 )  THEN
1804             pmf2a(:,1) = 0.0_wp
1805             pmf2b(:,1) = 0.0_wp
1806          ENDIF
1807          IF ( index_oc < 0 )  THEN
1808             pmf2a(:,2) = 0.0_wp
1809             pmf2b(:,2) = 0.0_wp
1810          ENDIF
1811          IF ( index_bc < 0 )  THEN
1812             pmf2a(:,3) = 0.0_wp
1813             pmf2b(:,3) = 0.0_wp
1814          ENDIF
1815          IF ( index_du < 0 )  THEN
1816             pmf2a(:,4) = 0.0_wp
1817             pmf2b(:,4) = 0.0_wp
1818          ENDIF
1819          IF ( index_ss < 0 )  THEN
1820             pmf2a(:,5) = 0.0_wp
1821             pmf2b(:,5) = 0.0_wp
1822          ENDIF
1823          IF ( index_no < 0 )  THEN
1824             pmf2a(:,6) = 0.0_wp
1825             pmf2b(:,6) = 0.0_wp
1826          ENDIF
1827          IF ( index_nh < 0 )  THEN
1828             pmf2a(:,7) = 0.0_wp
1829             pmf2b(:,7) = 0.0_wp
1830          ENDIF
1831
1832          IF ( SUM( pmf2a ) < 0.00001_wp  .AND.  SUM( pmf2b ) < 0.00001_wp )  THEN
1833             message_string = 'Error in initialising mass fractions of chemical components. ' //   &
1834                              'Check that all chemical components are included in parameter file!'
1835             CALL message( 'salsa_mod: aerosol_init', 'PA0606', 2, 2, 0, 6, 0 ) 
1836          ENDIF
1837!
1838!--       Then normalise the mass fraction so that SUM = 1
1839          DO  k = nzb, nzt+1
1840             pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) )
1841             IF ( SUM( pmf2b(k,:) ) > 0.0_wp )  pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) )
1842          ENDDO
1843
1844          DEALLOCATE( pr_z, pr_mass_fracs_a, pr_mass_fracs_b )
1845
1846       ELSE
1847          message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) //  &
1848                           ' for SALSA missing!'
1849          CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 )
1850!
1851!--       Close input file
1852          CALL close_input_file( id_dyn )
1853       ENDIF   ! netcdf_extend
1854
1855#else
1856       message_string = 'init_aerosol_type = 1 but preprocessor directive __netcdf is not used '// &
1857                        'in compiling!'
1858       CALL message( 'salsa_mod: aerosol_init', 'PA0608', 1, 2, 0, 6, 0 )
1859
1860#endif
1861
1862    ELSEIF ( init_aerosol_type == 0 )  THEN
1863!
1864!--    Mass fractions for species in a and b-bins
1865       IF ( index_so4 > 0 )  THEN
1866          pmf2a(:,1) = mass_fracs_a(index_so4)
1867          pmf2b(:,1) = mass_fracs_b(index_so4)
1868       ENDIF
1869       IF ( index_oc > 0 )  THEN
1870          pmf2a(:,2) = mass_fracs_a(index_oc)
1871          pmf2b(:,2) = mass_fracs_b(index_oc)
1872       ENDIF
1873       IF ( index_bc > 0 )  THEN
1874          pmf2a(:,3) = mass_fracs_a(index_bc)
1875          pmf2b(:,3) = mass_fracs_b(index_bc)
1876       ENDIF
1877       IF ( index_du > 0 )  THEN
1878          pmf2a(:,4) = mass_fracs_a(index_du)
1879          pmf2b(:,4) = mass_fracs_b(index_du)
1880       ENDIF
1881       IF ( index_ss > 0 )  THEN
1882          pmf2a(:,5) = mass_fracs_a(index_ss)
1883          pmf2b(:,5) = mass_fracs_b(index_ss)
1884       ENDIF
1885       IF ( index_no > 0 )  THEN
1886          pmf2a(:,6) = mass_fracs_a(index_no)
1887          pmf2b(:,6) = mass_fracs_b(index_no)
1888       ENDIF
1889       IF ( index_nh > 0 )  THEN
1890          pmf2a(:,7) = mass_fracs_a(index_nh)
1891          pmf2b(:,7) = mass_fracs_b(index_nh)
1892       ENDIF
1893       DO  k = nzb, nzt+1
1894          pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) )
1895          IF ( SUM( pmf2b(k,:) ) > 0.0_wp ) pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) )
1896       ENDDO
1897
1898       CALL size_distribution( n_lognorm, dpg, sigmag, nsect )
1899!
1900!--    Normalize by the given total number concentration
1901       nsect = nsect * SUM( n_lognorm ) / SUM( nsect )
1902       DO  ib = start_subrange_1a, end_subrange_2b
1903          pndist(:,ib) = nsect(ib)
1904       ENDDO
1905    ENDIF
1906
1907    IF ( init_gases_type == 1 )  THEN
1908!
1909!--    Read input profiles from PIDS_CHEM
1910#if defined( __netcdf )
1911!
1912!--    Location-dependent size distributions and compositions.
1913       INQUIRE( FILE = TRIM( input_file_dynamic ) //  TRIM( coupling_char ), EXIST = netcdf_extend )
1914       IF ( netcdf_extend  .AND.  .NOT. salsa_gases_from_chem )  THEN
1915!
1916!--       Open file in read-only mode
1917          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
1918!
1919!--       Inquire dimensions:
1920          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' )
1921          IF ( pr_nz /= nz )  THEN
1922             WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//&
1923                                        'the number of numeric grid points.'
1924             CALL message( 'aerosol_init', 'PA0609', 1, 2, 0, 6, 0 )
1925          ENDIF
1926!
1927!--       Read vertical profiles of gases:
1928          CALL get_variable( id_dyn, 'init_atmosphere_h2so4', salsa_gas(1)%init(nzb+1:nzt) )
1929          CALL get_variable( id_dyn, 'init_atmosphere_hno3',  salsa_gas(2)%init(nzb+1:nzt) )
1930          CALL get_variable( id_dyn, 'init_atmosphere_nh3',   salsa_gas(3)%init(nzb+1:nzt) )
1931          CALL get_variable( id_dyn, 'init_atmosphere_ocnv',  salsa_gas(4)%init(nzb+1:nzt) )
1932          CALL get_variable( id_dyn, 'init_atmosphere_ocsv',  salsa_gas(5)%init(nzb+1:nzt) )
1933!
1934!--       Set Neumann top and surface boundary condition for initial + initialise concentrations
1935          DO  ig = 1, ngases_salsa
1936             salsa_gas(ig)%init(nzb)   =  salsa_gas(ig)%init(nzb+1)
1937             salsa_gas(ig)%init(nzt+1) =  salsa_gas(ig)%init(nzt)
1938             IF ( .NOT. read_restart_data_salsa )  THEN
1939                DO  k = nzb, nzt+1
1940                   salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
1941                ENDDO
1942             ENDIF
1943          ENDDO
1944
1945       ELSEIF ( .NOT. netcdf_extend  .AND.  .NOT.  salsa_gases_from_chem )  THEN
1946          message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) //  &
1947                           ' for SALSA missing!'
1948          CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 )
1949!
1950!--       Close input file
1951          CALL close_input_file( id_dyn )
1952       ENDIF   ! netcdf_extend
1953#else
1954       message_string = 'init_gases_type = 1 but preprocessor directive __netcdf is not used in '//&
1955                        'compiling!'
1956       CALL message( 'salsa_mod: aerosol_init', 'PA0611', 1, 2, 0, 6, 0 )
1957
1958#endif
1959
1960    ENDIF
1961!
1962!-- Both SO4 and OC are included, so use the given mass fractions
1963    IF ( index_oc > 0  .AND.  index_so4 > 0 )  THEN
1964       pmfoc1a(:) = pmf2a(:,2) / ( pmf2a(:,2) + pmf2a(:,1) )  ! Normalize
1965!
1966!-- Pure organic carbon
1967    ELSEIF ( index_oc > 0 )  THEN
1968       pmfoc1a(:) = 1.0_wp
1969!
1970!-- Pure SO4
1971    ELSEIF ( index_so4 > 0 )  THEN
1972       pmfoc1a(:) = 0.0_wp
1973
1974    ELSE
1975       message_string = 'Either OC or SO4 must be active for aerosol region 1a!'
1976       CALL message( 'salsa_mod: aerosol_init', 'PA0612', 1, 2, 0, 6, 0 )
1977    ENDIF
1978
1979!
1980!-- Initialize concentrations
1981    DO  i = nxlg, nxrg
1982       DO  j = nysg, nyng
1983          DO  k = nzb, nzt+1
1984!
1985!--          Predetermine flag to mask topography
1986             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1987!
1988!--          a) Number concentrations
1989!--          Region 1:
1990             DO  ib = start_subrange_1a, end_subrange_1a
1991                IF ( .NOT. read_restart_data_salsa )  THEN
1992                   aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
1993                ENDIF
1994                IF ( prunmode == 1 )  THEN
1995                   aerosol_number(ib)%init = pndist(:,ib)
1996                ENDIF
1997             ENDDO
1998!
1999!--          Region 2:
2000             IF ( nreg > 1 )  THEN
2001                DO  ib = start_subrange_2a, end_subrange_2a
2002                   IF ( .NOT. read_restart_data_salsa )  THEN
2003                      aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
2004                   ENDIF
2005                   IF ( prunmode == 1 )  THEN
2006                      aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib)
2007                   ENDIF
2008                ENDDO
2009                IF ( .NOT. no_insoluble )  THEN
2010                   DO  ib = start_subrange_2b, end_subrange_2b
2011                      IF ( pnf2a(k) < 1.0_wp )  THEN
2012                         IF ( .NOT. read_restart_data_salsa )  THEN
2013                            aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *    &
2014                                                             pndist(k,ib) * flag
2015                         ENDIF
2016                         IF ( prunmode == 1 )  THEN
2017                            aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp - nf2a ) * pndist(:,ib)
2018                         ENDIF
2019                      ENDIF
2020                   ENDDO
2021                ENDIF
2022             ENDIF
2023!
2024!--          b) Aerosol mass concentrations
2025!--             bin subrange 1: done here separately due to the SO4/OC convention
2026!
2027!--          SO4:
2028             IF ( index_so4 > 0 )  THEN
2029                ss = ( index_so4 - 1 ) * nbins_aerosol + start_subrange_1a !< start
2030                ee = ( index_so4 - 1 ) * nbins_aerosol + end_subrange_1a !< end
2031                ib = start_subrange_1a
2032                DO  ic = ss, ee
2033                   IF ( .NOT. read_restart_data_salsa )  THEN
2034                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) *          &
2035                                                     pndist(k,ib) * core(ib) * arhoh2so4 * flag
2036                   ENDIF
2037                   IF ( prunmode == 1 )  THEN
2038                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib) &
2039                                                 * core(ib) * arhoh2so4
2040                   ENDIF
2041                   ib = ib+1
2042                ENDDO
2043             ENDIF
2044!
2045!--          OC:
2046             IF ( index_oc > 0 ) THEN
2047                ss = ( index_oc - 1 ) * nbins_aerosol + start_subrange_1a !< start
2048                ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a !< end
2049                ib = start_subrange_1a
2050                DO  ic = ss, ee
2051                   IF ( .NOT. read_restart_data_salsa )  THEN
2052                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *    &
2053                                                     core(ib) * arhooc * flag
2054                   ENDIF
2055                   IF ( prunmode == 1 )  THEN
2056                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *        &
2057                                                 core(ib) * arhooc
2058                   ENDIF
2059                   ib = ib+1
2060                ENDDO
2061             ENDIF
2062          ENDDO !< k
2063
2064          prunmode = 3  ! Init only once
2065
2066       ENDDO !< j
2067    ENDDO !< i
2068
2069!
2070!-- c) Aerosol mass concentrations
2071!--    bin subrange 2:
2072    IF ( nreg > 1 ) THEN
2073
2074       IF ( index_so4 > 0 ) THEN
2075          CALL set_aero_mass( index_so4, pmf2a(:,1), pmf2b(:,1), pnf2a, pndist, core, arhoh2so4 )
2076       ENDIF
2077       IF ( index_oc > 0 ) THEN
2078          CALL set_aero_mass( index_oc, pmf2a(:,2), pmf2b(:,2), pnf2a, pndist, core, arhooc )
2079       ENDIF
2080       IF ( index_bc > 0 ) THEN
2081          CALL set_aero_mass( index_bc, pmf2a(:,3), pmf2b(:,3), pnf2a, pndist, core, arhobc )
2082       ENDIF
2083       IF ( index_du > 0 ) THEN
2084          CALL set_aero_mass( index_du, pmf2a(:,4), pmf2b(:,4), pnf2a, pndist, core, arhodu )
2085       ENDIF
2086       IF ( index_ss > 0 ) THEN
2087          CALL set_aero_mass( index_ss, pmf2a(:,5), pmf2b(:,5), pnf2a, pndist, core, arhoss )
2088       ENDIF
2089       IF ( index_no > 0 ) THEN
2090          CALL set_aero_mass( index_no, pmf2a(:,6), pmf2b(:,6), pnf2a, pndist, core, arhohno3 )
2091       ENDIF
2092       IF ( index_nh > 0 ) THEN
2093          CALL set_aero_mass( index_nh, pmf2a(:,7), pmf2b(:,7), pnf2a, pndist, core, arhonh3 )
2094       ENDIF
2095
2096    ENDIF
2097
2098 END SUBROUTINE aerosol_init
2099
2100!------------------------------------------------------------------------------!
2101! Description:
2102! ------------
2103!> Create a lognormal size distribution and discretise to a sectional
2104!> representation.
2105!------------------------------------------------------------------------------!
2106 SUBROUTINE size_distribution( in_ntot, in_dpg, in_sigma, psd_sect )
2107
2108    IMPLICIT NONE
2109
2110    INTEGER(iwp) ::  ib         !< running index: bin
2111    INTEGER(iwp) ::  iteration  !< running index: iteration
2112
2113    REAL(wp) ::  d1         !< particle diameter (m, dummy)
2114    REAL(wp) ::  d2         !< particle diameter (m, dummy)
2115    REAL(wp) ::  delta_d    !< (d2-d1)/10
2116    REAL(wp) ::  deltadp    !< bin width
2117    REAL(wp) ::  dmidi      !< ( d1 + d2 ) / 2
2118
2119    REAL(wp), DIMENSION(:), INTENT(in) ::  in_dpg    !< geometric mean diameter (m)
2120    REAL(wp), DIMENSION(:), INTENT(in) ::  in_ntot   !< number conc. (#/m3)
2121    REAL(wp), DIMENSION(:), INTENT(in) ::  in_sigma  !< standard deviation
2122
2123    REAL(wp), DIMENSION(:), INTENT(inout) ::  psd_sect  !< sectional size distribution
2124
2125    DO  ib = start_subrange_1a, end_subrange_2b
2126       psd_sect(ib) = 0.0_wp
2127!
2128!--    Particle diameter at the low limit (largest in the bin) (m)
2129       d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp
2130!
2131!--    Particle diameter at the high limit (smallest in the bin) (m)
2132       d2 = ( aero(ib)%vhilim / api6 )**0.33333333_wp
2133!
2134!--    Span of particle diameter in a bin (m)
2135       delta_d = 0.1_wp * ( d2 - d1 )
2136!
2137!--    Iterate:
2138       DO  iteration = 1, 10
2139          d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp + ( ib - 1) * delta_d
2140          d2 = d1 + delta_d
2141          dmidi = 0.5_wp * ( d1 + d2 )
2142          deltadp = LOG10( d2 / d1 )
2143!
2144!--       Size distribution
2145!--       in_ntot = total number, total area, or total volume concentration
2146!--       in_dpg = geometric-mean number, area, or volume diameter
2147!--       n(k) = number, area, or volume concentration in a bin
2148          psd_sect(ib) = psd_sect(ib) + SUM( in_ntot * deltadp / ( SQRT( 2.0_wp * pi ) *           &
2149                        LOG10( in_sigma ) ) * EXP( -LOG10( dmidi / in_dpg )**2.0_wp /              &
2150                        ( 2.0_wp * LOG10( in_sigma ) ** 2.0_wp ) ) )
2151
2152       ENDDO
2153    ENDDO
2154
2155 END SUBROUTINE size_distribution
2156
2157!------------------------------------------------------------------------------!
2158! Description:
2159! ------------
2160!> Sets the mass concentrations to aerosol arrays in 2a and 2b.
2161!>
2162!> Tomi Raatikainen, FMI, 29.2.2016
2163!------------------------------------------------------------------------------!
2164 SUBROUTINE set_aero_mass( ispec, pmf2a, pmf2b, pnf2a, pndist, pcore, prho )
2165
2166    IMPLICIT NONE
2167
2168    INTEGER(iwp) ::  ee        !< index: end
2169    INTEGER(iwp) ::  i         !< loop index
2170    INTEGER(iwp) ::  ib        !< loop index
2171    INTEGER(iwp) ::  ic        !< loop index
2172    INTEGER(iwp) ::  j         !< loop index
2173    INTEGER(iwp) ::  k         !< loop index
2174    INTEGER(iwp) ::  prunmode  !< 1 = initialise
2175    INTEGER(iwp) ::  ss        !< index: start
2176
2177    INTEGER(iwp), INTENT(in) :: ispec  !< Aerosol species index
2178
2179    REAL(wp) ::  flag   !< flag to mask topography grid points
2180
2181    REAL(wp), INTENT(in) ::  prho !< Aerosol density
2182
2183    REAL(wp), DIMENSION(nbins_aerosol), INTENT(in) ::  pcore !< Aerosol bin mid core volume
2184    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pnf2a !< Number fraction for 2a
2185    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pmf2a !< Mass distributions for a
2186    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pmf2b !< and b bins
2187
2188    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol), INTENT(in) ::  pndist !< Aerosol size distribution
2189
2190    prunmode = 1
2191
2192    DO i = nxlg, nxrg
2193       DO j = nysg, nyng
2194          DO k = nzb, nzt+1
2195!
2196!--          Predetermine flag to mask topography
2197             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
2198!
2199!--          Regime 2a:
2200             ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2a
2201             ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2a
2202             ib = start_subrange_2a
2203             DO ic = ss, ee
2204                IF ( .NOT. read_restart_data_salsa )  THEN
2205                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib)&
2206                                                  * pcore(ib) * prho * flag
2207                ENDIF
2208                IF ( prunmode == 1 )  THEN
2209                   aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) *  &
2210                                              pcore(ib) * prho
2211                ENDIF
2212                ib = ib + 1
2213             ENDDO
2214!
2215!--          Regime 2b:
2216             IF ( .NOT. no_insoluble )  THEN
2217                ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2b
2218                ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2b
2219                ib = start_subrange_2a
2220                DO ic = ss, ee
2221                   IF ( .NOT. read_restart_data_salsa )  THEN
2222                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k))&
2223                                                     * pndist(k,ib) * pcore(ib) * prho * flag
2224                   ENDIF
2225                   IF ( prunmode == 1 )  THEN
2226                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) * &
2227                                                 pndist(k,ib) * pcore(ib) * prho
2228                   ENDIF
2229                   ib = ib + 1
2230                ENDDO  ! c
2231
2232             ENDIF
2233          ENDDO   ! k
2234
2235          prunmode = 3  ! Init only once
2236
2237       ENDDO   ! j
2238    ENDDO   ! i
2239
2240 END SUBROUTINE set_aero_mass
2241
2242!------------------------------------------------------------------------------!
2243! Description:
2244! ------------
2245!> Initialise the matching between surface types in LSM and deposition models.
2246!> Do the matching based on Zhang et al. (2001). Atmos. Environ. 35, 549-560
2247!> (here referred as Z01).
2248!------------------------------------------------------------------------------!
2249 SUBROUTINE init_deposition
2250
2251    USE surface_mod,                                                                               &
2252        ONLY:  surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
2253
2254    IMPLICIT NONE
2255
2256    INTEGER(iwp) ::  l  !< loop index for vertical surfaces
2257
2258    LOGICAL :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM surfaces)
2259
2260    IF ( depo_pcm_par == 'zhang2001' )  THEN
2261       depo_pcm_par_num = 1
2262    ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
2263       depo_pcm_par_num = 2
2264    ENDIF
2265
2266    IF ( depo_surf_par == 'zhang2001' )  THEN
2267       depo_surf_par_num = 1
2268    ELSEIF ( depo_surf_par == 'petroff2010' )  THEN
2269       depo_surf_par_num = 2
2270    ENDIF
2271!
2272!-- LSM: Pavement, vegetation and water
2273    IF ( nldepo_surf  .AND.  land_surface )  THEN
2274       match_lsm = .TRUE.
2275       ALLOCATE( lsm_to_depo_h%match_lupg(1:surf_lsm_h%ns),                                         &
2276                 lsm_to_depo_h%match_luvw(1:surf_lsm_h%ns),                                         &
2277                 lsm_to_depo_h%match_luww(1:surf_lsm_h%ns) )
2278       lsm_to_depo_h%match_lupg = 0
2279       lsm_to_depo_h%match_luvw = 0
2280       lsm_to_depo_h%match_luww = 0
2281       CALL match_sm_zhang( surf_lsm_h, lsm_to_depo_h%match_lupg, lsm_to_depo_h%match_luvw,        &
2282                            lsm_to_depo_h%match_luww, match_lsm )
2283       DO  l = 0, 3
2284          ALLOCATE( lsm_to_depo_v(l)%match_lupg(1:surf_lsm_v(l)%ns),                               &
2285                    lsm_to_depo_v(l)%match_luvw(1:surf_lsm_v(l)%ns),                               &
2286                    lsm_to_depo_v(l)%match_luww(1:surf_lsm_v(l)%ns) )
2287          lsm_to_depo_v(l)%match_lupg = 0
2288          lsm_to_depo_v(l)%match_luvw = 0
2289          lsm_to_depo_v(l)%match_luww = 0
2290          CALL match_sm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match_lupg,                         &
2291                               lsm_to_depo_v(l)%match_luvw, lsm_to_depo_v(l)%match_luww, match_lsm )
2292       ENDDO
2293    ENDIF
2294!
2295!-- USM: Green roofs/walls, wall surfaces and windows
2296    IF ( nldepo_surf  .AND.  urban_surface )  THEN
2297       match_lsm = .FALSE.
2298       ALLOCATE( usm_to_depo_h%match_lupg(1:surf_usm_h%ns),                                        &
2299                 usm_to_depo_h%match_luvw(1:surf_usm_h%ns),                                        &
2300                 usm_to_depo_h%match_luww(1:surf_usm_h%ns) )
2301       usm_to_depo_h%match_lupg = 0
2302       usm_to_depo_h%match_luvw = 0
2303       usm_to_depo_h%match_luww = 0
2304       CALL match_sm_zhang( surf_usm_h, usm_to_depo_h%match_lupg, usm_to_depo_h%match_luvw,        &
2305                            usm_to_depo_h%match_luww, match_lsm )
2306       DO  l = 0, 3
2307          ALLOCATE( usm_to_depo_v(l)%match_lupg(1:surf_usm_v(l)%ns),                               &
2308                    usm_to_depo_v(l)%match_luvw(1:surf_usm_v(l)%ns),                               &
2309                    usm_to_depo_v(l)%match_luww(1:surf_usm_v(l)%ns) )
2310          usm_to_depo_v(l)%match_lupg = 0
2311          usm_to_depo_v(l)%match_luvw = 0
2312          usm_to_depo_v(l)%match_luww = 0
2313          CALL match_sm_zhang( surf_usm_v(l), usm_to_depo_v(l)%match_lupg,                         &
2314                               usm_to_depo_v(l)%match_luvw, usm_to_depo_v(l)%match_luww, match_lsm )
2315       ENDDO
2316    ENDIF
2317
2318    IF ( nldepo_pcm )  THEN
2319       SELECT CASE ( depo_pcm_type )
2320          CASE ( 'evergreen_needleleaf' )
2321             depo_pcm_type_num = 1
2322          CASE ( 'evergreen_broadleaf' )
2323             depo_pcm_type_num = 2
2324          CASE ( 'deciduous_needleleaf' )
2325             depo_pcm_type_num = 3
2326          CASE ( 'deciduous_broadleaf' )
2327             depo_pcm_type_num = 4
2328          CASE DEFAULT
2329             message_string = 'depo_pcm_type not set correctly.'
2330             CALL message( 'salsa_mod: init_deposition', 'PA0613', 1, 2, 0, 6, 0 )
2331       END SELECT
2332    ENDIF
2333
2334 END SUBROUTINE init_deposition
2335
2336!------------------------------------------------------------------------------!
2337! Description:
2338! ------------
2339!> Match the surface types in PALM and Zhang et al. 2001 deposition module
2340!------------------------------------------------------------------------------!
2341 SUBROUTINE match_sm_zhang( surf, match_pav_green, match_veg_wall, match_wat_win, match_lsm )
2342
2343    USE surface_mod,                                                           &
2344        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
2345
2346    IMPLICIT NONE
2347
2348    INTEGER(iwp) ::  m              !< index for surface elements
2349    INTEGER(iwp) ::  pav_type_palm  !< pavement / green wall type in PALM
2350    INTEGER(iwp) ::  veg_type_palm  !< vegetation / wall type in PALM
2351    INTEGER(iwp) ::  wat_type_palm  !< water / window type in PALM
2352
2353    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_pav_green  !<  matching pavement/green walls
2354    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_veg_wall   !<  matching vegetation/walls
2355    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_wat_win    !<  matching water/windows
2356
2357    LOGICAL, INTENT(in) :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM)
2358
2359    TYPE(surf_type), INTENT(in) :: surf  !< respective surface type
2360
2361    DO  m = 1, surf%ns
2362       IF ( match_lsm )  THEN
2363!
2364!--       Vegetation (LSM):
2365          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
2366             veg_type_palm = surf%vegetation_type(m)
2367             SELECT CASE ( veg_type_palm )
2368                CASE ( 0 )
2369                   message_string = 'No vegetation type defined.'
2370                   CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
2371                CASE ( 1 )  ! bare soil
2372                   match_veg_wall(m) = 6  ! grass in Z01
2373                CASE ( 2 )  ! crops, mixed farming
2374                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
2375                CASE ( 3 )  ! short grass
2376                   match_veg_wall(m) = 6  ! grass in Z01
2377                CASE ( 4 )  ! evergreen needleleaf trees
2378                    match_veg_wall(m) = 1  ! evergreen needleleaf trees in Z01
2379                CASE ( 5 )  ! deciduous needleleaf trees
2380                   match_veg_wall(m) = 3  ! deciduous needleleaf trees in Z01
2381                CASE ( 6 )  ! evergreen broadleaf trees
2382                   match_veg_wall(m) = 2  ! evergreen broadleaf trees in Z01
2383                CASE ( 7 )  ! deciduous broadleaf trees
2384                   match_veg_wall(m) = 4  ! deciduous broadleaf trees in Z01
2385                CASE ( 8 )  ! tall grass
2386                   match_veg_wall(m) = 6  ! grass in Z01
2387                CASE ( 9 )  ! desert
2388                   match_veg_wall(m) = 8  ! desert in Z01
2389                CASE ( 10 )  ! tundra
2390                   match_veg_wall(m) = 9  ! tundra in Z01
2391                CASE ( 11 )  ! irrigated crops
2392                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
2393                CASE ( 12 )  ! semidesert
2394                   match_veg_wall(m) = 8  ! desert in Z01
2395                CASE ( 13 )  ! ice caps and glaciers
2396                   match_veg_wall(m) = 12  ! ice cap and glacier in Z01
2397                CASE ( 14 )  ! bogs and marshes
2398                   match_veg_wall(m) = 11  ! wetland with plants in Z01
2399                CASE ( 15 )  ! evergreen shrubs
2400                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2401                CASE ( 16 )  ! deciduous shrubs
2402                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2403                CASE ( 17 )  ! mixed forest/woodland
2404                   match_veg_wall(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
2405                CASE ( 18 )  ! interrupted forest
2406                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2407             END SELECT
2408          ENDIF
2409!
2410!--       Pavement (LSM):
2411          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
2412             pav_type_palm = surf%pavement_type(m)
2413             IF ( pav_type_palm == 0 )  THEN  ! error
2414                message_string = 'No pavement type defined.'
2415                CALL message( 'salsa_mod: match_sm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
2416             ELSE
2417                match_pav_green(m) = 15  ! urban in Z01
2418             ENDIF
2419          ENDIF
2420!
2421!--       Water (LSM):
2422          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
2423             wat_type_palm = surf%water_type(m)
2424             IF ( wat_type_palm == 0 )  THEN  ! error
2425                message_string = 'No water type defined.'
2426                CALL message( 'salsa_mod: match_sm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
2427             ELSEIF ( wat_type_palm == 3 )  THEN
2428                match_wat_win(m) = 14  ! ocean in Z01
2429             ELSEIF ( wat_type_palm == 1  .OR.  wat_type_palm == 2 .OR.  wat_type_palm == 4        &
2430                      .OR.  wat_type_palm == 5  )  THEN
2431                match_wat_win(m) = 13  ! inland water in Z01
2432             ENDIF
2433          ENDIF
2434       ELSE
2435!
2436!--       Wall surfaces (USM):
2437          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
2438             match_veg_wall(m) = 15  ! urban in Z01
2439          ENDIF
2440!
2441!--       Green walls and roofs (USM):
2442          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
2443             match_pav_green(m) =  6 ! (short) grass in Z01
2444          ENDIF
2445!
2446!--       Windows (USM):
2447          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
2448             match_wat_win(m) = 15  ! urban in Z01
2449          ENDIF
2450       ENDIF
2451
2452    ENDDO
2453
2454 END SUBROUTINE match_sm_zhang
2455
2456!------------------------------------------------------------------------------!
2457! Description:
2458! ------------
2459!> Swapping of timelevels
2460!------------------------------------------------------------------------------!
2461 SUBROUTINE salsa_swap_timelevel( mod_count )
2462
2463    IMPLICIT NONE
2464
2465    INTEGER(iwp) ::  ib   !<
2466    INTEGER(iwp) ::  ic   !<
2467    INTEGER(iwp) ::  icc  !<
2468    INTEGER(iwp) ::  ig   !<
2469
2470    INTEGER(iwp), INTENT(IN) ::  mod_count  !<
2471
2472    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
2473
2474       SELECT CASE ( mod_count )
2475
2476          CASE ( 0 )
2477
2478             DO  ib = 1, nbins_aerosol
2479                aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => nconc_1(:,:,:,ib)
2480                aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,ib)
2481
2482                DO  ic = 1, ncomponents_mass
2483                   icc = ( ic-1 ) * nbins_aerosol + ib
2484                   aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => mconc_1(:,:,:,icc)
2485                   aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,icc)
2486                ENDDO
2487             ENDDO
2488
2489             IF ( .NOT. salsa_gases_from_chem )  THEN
2490                DO  ig = 1, ngases_salsa
2491                   salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => gconc_1(:,:,:,ig)
2492                   salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,ig)
2493                ENDDO
2494             ENDIF
2495
2496          CASE ( 1 )
2497
2498             DO  ib = 1, nbins_aerosol
2499                aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => nconc_2(:,:,:,ib)
2500                aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,ib)
2501                DO  ic = 1, ncomponents_mass
2502                   icc = ( ic-1 ) * nbins_aerosol + ib
2503                   aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => mconc_2(:,:,:,icc)
2504                   aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,icc)
2505                ENDDO
2506             ENDDO
2507
2508             IF ( .NOT. salsa_gases_from_chem )  THEN
2509                DO  ig = 1, ngases_salsa
2510                   salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => gconc_2(:,:,:,ig)
2511                   salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,ig)
2512                ENDDO
2513             ENDIF
2514
2515       END SELECT
2516
2517    ENDIF
2518
2519 END SUBROUTINE salsa_swap_timelevel
2520
2521
2522!------------------------------------------------------------------------------!
2523! Description:
2524! ------------
2525!> This routine reads the respective restart data.
2526!------------------------------------------------------------------------------!
2527 SUBROUTINE salsa_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,      &
2528                             nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
2529
2530    IMPLICIT NONE
2531
2532    INTEGER(iwp) ::  ib              !<
2533    INTEGER(iwp) ::  ic              !<
2534    INTEGER(iwp) ::  ig              !<
2535    INTEGER(iwp) ::  k               !<
2536    INTEGER(iwp) ::  nxlc            !<
2537    INTEGER(iwp) ::  nxlf            !<
2538    INTEGER(iwp) ::  nxl_on_file     !<
2539    INTEGER(iwp) ::  nxrc            !<
2540    INTEGER(iwp) ::  nxrf            !<
2541    INTEGER(iwp) ::  nxr_on_file     !<
2542    INTEGER(iwp) ::  nync            !<
2543    INTEGER(iwp) ::  nynf            !<
2544    INTEGER(iwp) ::  nyn_on_file     !<
2545    INTEGER(iwp) ::  nysc            !<
2546    INTEGER(iwp) ::  nysf            !<
2547    INTEGER(iwp) ::  nys_on_file     !<
2548
2549    LOGICAL, INTENT(OUT)  ::  found  !<
2550
2551    REAL(wp), &
2552       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
2553
2554    found = .FALSE.
2555
2556    IF ( read_restart_data_salsa )  THEN
2557
2558       SELECT CASE ( restart_string(1:length) )
2559
2560          CASE ( 'aerosol_number' )
2561             DO  ib = 1, nbins_aerosol
2562                IF ( k == 1 )  READ ( 13 ) tmp_3d
2563                aerosol_number(ib)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
2564                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2565                found = .TRUE.
2566             ENDDO
2567
2568          CASE ( 'aerosol_mass' )
2569             DO  ic = 1, ncomponents_mass * nbins_aerosol
2570                IF ( k == 1 )  READ ( 13 ) tmp_3d
2571                aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
2572                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2573                found = .TRUE.
2574             ENDDO
2575
2576          CASE ( 'salsa_gas' )
2577             DO  ig = 1, ngases_salsa
2578                IF ( k == 1 )  READ ( 13 ) tmp_3d
2579                salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                    &
2580                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2581                found = .TRUE.
2582             ENDDO
2583
2584          CASE DEFAULT
2585             found = .FALSE.
2586
2587       END SELECT
2588    ENDIF
2589
2590 END SUBROUTINE salsa_rrd_local
2591
2592!------------------------------------------------------------------------------!
2593! Description:
2594! ------------
2595!> This routine writes the respective restart data.
2596!> Note that the following input variables in PARIN have to be equal between
2597!> restart runs:
2598!>    listspec, nbin, nbin2, nf2a, ncc, mass_fracs_a, mass_fracs_b
2599!------------------------------------------------------------------------------!
2600 SUBROUTINE salsa_wrd_local
2601
2602    IMPLICIT NONE
2603
2604    INTEGER(iwp) ::  ib   !<
2605    INTEGER(iwp) ::  ic   !<
2606    INTEGER(iwp) ::  ig  !<
2607
2608    IF ( write_binary  .AND.  write_binary_salsa )  THEN
2609
2610       CALL wrd_write_string( 'aerosol_number' )
2611       DO  ib = 1, nbins_aerosol
2612          WRITE ( 14 )  aerosol_number(ib)%conc
2613       ENDDO
2614
2615       CALL wrd_write_string( 'aerosol_mass' )
2616       DO  ic = 1, nbins_aerosol * ncomponents_mass
2617          WRITE ( 14 )  aerosol_mass(ic)%conc
2618       ENDDO
2619
2620       CALL wrd_write_string( 'salsa_gas' )
2621       DO  ig = 1, ngases_salsa
2622          WRITE ( 14 )  salsa_gas(ig)%conc
2623       ENDDO
2624
2625    ENDIF
2626
2627 END SUBROUTINE salsa_wrd_local
2628
2629!------------------------------------------------------------------------------!
2630! Description:
2631! ------------
2632!> Performs necessary unit and dimension conversion between the host model and
2633!> SALSA module, and calls the main SALSA routine.
2634!> Partially adobted form the original SALSA boxmodel version.
2635!> Now takes masses in as kg/kg from LES!! Converted to m3/m3 for SALSA
2636!> 05/2016 Juha: This routine is still pretty much in its original shape.
2637!>               It's dumb as a mule and twice as ugly, so implementation of
2638!>               an improved solution is necessary sooner or later.
2639!> Juha Tonttila, FMI, 2014
2640!> Jaakko Ahola, FMI, 2016
2641!> Only aerosol processes included, Mona Kurppa, UHel, 2017
2642!------------------------------------------------------------------------------!
2643 SUBROUTINE salsa_driver( i, j, prunmode )
2644
2645    USE arrays_3d,                                                                                 &
2646        ONLY: pt_p, q_p, u, v, w
2647
2648    USE plant_canopy_model_mod,                                                                    &
2649        ONLY: lad_s
2650
2651    USE surface_mod,                                                                               &
2652        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
2653
2654    IMPLICIT NONE
2655
2656    INTEGER(iwp) ::  endi    !< end index
2657    INTEGER(iwp) ::  ib      !< loop index
2658    INTEGER(iwp) ::  ic      !< loop index
2659    INTEGER(iwp) ::  ig      !< loop index
2660    INTEGER(iwp) ::  k_wall  !< vertical index of topography top
2661    INTEGER(iwp) ::  k       !< loop index
2662    INTEGER(iwp) ::  l       !< loop index
2663    INTEGER(iwp) ::  nc_h2o  !< index of H2O in the prtcl index table
2664    INTEGER(iwp) ::  ss      !< loop index
2665    INTEGER(iwp) ::  str     !< start index
2666    INTEGER(iwp) ::  vc      !< default index in prtcl
2667
2668    INTEGER(iwp), INTENT(in) ::  i         !< loop index
2669    INTEGER(iwp), INTENT(in) ::  j         !< loop index
2670    INTEGER(iwp), INTENT(in) ::  prunmode  !< 1: Initialization, 2: Spinup, 3: Regular runtime
2671
2672    REAL(wp) ::  cw_old  !< previous H2O mixing ratio
2673    REAL(wp) ::  flag    !< flag to mask topography grid points
2674    REAL(wp) ::  in_lad  !< leaf area density (m2/m3)
2675    REAL(wp) ::  in_rh   !< relative humidity
2676    REAL(wp) ::  zgso4   !< SO4
2677    REAL(wp) ::  zghno3  !< HNO3
2678    REAL(wp) ::  zgnh3   !< NH3
2679    REAL(wp) ::  zgocnv  !< non-volatile OC
2680    REAL(wp) ::  zgocsv  !< semi-volatile OC
2681
2682    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_adn  !< air density (kg/m3)
2683    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cs   !< H2O sat. vapour conc.
2684    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cw   !< H2O vapour concentration
2685    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_p    !< pressure (Pa)
2686    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_t    !< temperature (K)
2687    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_u    !< wind magnitude (m/s)
2688    REAL(wp), DIMENSION(nzb:nzt+1) ::  kvis    !< kinematic viscosity of air(m2/s)
2689    REAL(wp), DIMENSION(nzb:nzt+1) ::  ppm_to_nconc  !< Conversion factor from ppm to #/m3
2690
2691    REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) ::  schmidt_num  !< particle Schmidt number
2692    REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) ::  vd           !< particle fall seed (m/s)
2693
2694    TYPE(t_section), DIMENSION(nbins_aerosol) ::  lo_aero   !< additional variable for OpenMP
2695    TYPE(t_section), DIMENSION(nbins_aerosol) ::  aero_old  !< helper array
2696
2697    aero_old(:)%numc = 0.0_wp
2698    in_lad           = 0.0_wp
2699    in_u             = 0.0_wp
2700    kvis             = 0.0_wp
2701    lo_aero          = aero
2702    schmidt_num      = 0.0_wp
2703    vd               = 0.0_wp
2704    zgso4            = nclim
2705    zghno3           = nclim
2706    zgnh3            = nclim
2707    zgocnv           = nclim
2708    zgocsv           = nclim
2709!
2710!-- Aerosol number is always set, but mass can be uninitialized
2711    DO ib = 1, nbins_aerosol
2712       lo_aero(ib)%volc(:)  = 0.0_wp
2713       aero_old(ib)%volc(:) = 0.0_wp
2714    ENDDO
2715!
2716!-- Set the salsa runtime config (How to make this more efficient?)
2717    CALL set_salsa_runtime( prunmode )
2718!
2719!-- Calculate thermodynamic quantities needed in SALSA
2720    CALL salsa_thrm_ij( i, j, p_ij=in_p, temp_ij=in_t, cw_ij=in_cw, cs_ij=in_cs, adn_ij=in_adn )
2721!
2722!-- Magnitude of wind: needed for deposition
2723    IF ( lsdepo )  THEN
2724       in_u(nzb+1:nzt) = SQRT( ( 0.5_wp * ( u(nzb+1:nzt,j,i) + u(nzb+1:nzt,j,i+1) ) )**2 +         &
2725                               ( 0.5_wp * ( v(nzb+1:nzt,j,i) + v(nzb+1:nzt,j+1,i) ) )**2 +         &
2726                               ( 0.5_wp * ( w(nzb:nzt-1,j,i) + w(nzb+1:nzt,j,  i) ) )**2 )
2727    ENDIF
2728!
2729!-- Calculate conversion factors for gas concentrations
2730    ppm_to_nconc(:) = for_ppm_to_nconc * in_p(:) / in_t(:)
2731!
2732!-- Determine topography-top index on scalar grid
2733    k_wall = k_topo_top(j,i)
2734
2735    DO k = nzb+1, nzt
2736!
2737!--    Predetermine flag to mask topography
2738       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2739!
2740!--    Wind velocity for dry depositon on vegetation
2741       IF ( lsdepo_pcm  .AND.  plant_canopy )  THEN
2742          in_lad = lad_s( MAX( k-k_wall,0 ),j,i)
2743       ENDIF
2744!
2745!--    For initialization and spinup, limit the RH with the parameter rhlim
2746       IF ( prunmode < 3 ) THEN
2747          in_cw(k) = MIN( in_cw(k), in_cs(k) * rhlim )
2748       ELSE
2749          in_cw(k) = in_cw(k)
2750       ENDIF
2751       cw_old = in_cw(k) !* in_adn(k)
2752!
2753!--    Set volume concentrations:
2754!--    Sulphate (SO4) or sulphuric acid H2SO4
2755       IF ( index_so4 > 0 )  THEN
2756          vc = 1
2757          str = ( index_so4-1 ) * nbins_aerosol + 1    ! start index
2758          endi = index_so4 * nbins_aerosol             ! end index
2759          ic = 1
2760          DO ss = str, endi
2761             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2so4
2762             ic = ic+1
2763          ENDDO
2764          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2765       ENDIF
2766!
2767!--    Organic carbon (OC) compounds
2768       IF ( index_oc > 0 )  THEN
2769          vc = 2
2770          str = ( index_oc-1 ) * nbins_aerosol + 1
2771          endi = index_oc * nbins_aerosol
2772          ic = 1
2773          DO ss = str, endi
2774             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhooc
2775             ic = ic+1
2776          ENDDO
2777          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2778       ENDIF
2779!
2780!--    Black carbon (BC)
2781       IF ( index_bc > 0 )  THEN
2782          vc = 3
2783          str = ( index_bc-1 ) * nbins_aerosol + 1 + end_subrange_1a
2784          endi = index_bc * nbins_aerosol
2785          ic = 1 + end_subrange_1a
2786          DO ss = str, endi
2787             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhobc
2788             ic = ic+1
2789          ENDDO
2790          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2791       ENDIF
2792!
2793!--    Dust (DU)
2794       IF ( index_du > 0 )  THEN
2795          vc = 4
2796          str = ( index_du-1 ) * nbins_aerosol + 1 + end_subrange_1a
2797          endi = index_du * nbins_aerosol
2798          ic = 1 + end_subrange_1a
2799          DO ss = str, endi
2800             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhodu
2801             ic = ic+1
2802          ENDDO
2803          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2804       ENDIF
2805!
2806!--    Sea salt (SS)
2807       IF ( index_ss > 0 )  THEN
2808          vc = 5
2809          str = ( index_ss-1 ) * nbins_aerosol + 1 + end_subrange_1a
2810          endi = index_ss * nbins_aerosol
2811          ic = 1 + end_subrange_1a
2812          DO ss = str, endi
2813             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoss
2814             ic = ic+1
2815          ENDDO
2816          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2817       ENDIF
2818!
2819!--    Nitrate (NO(3-)) or nitric acid HNO3
2820       IF ( index_no > 0 )  THEN
2821          vc = 6
2822          str = ( index_no-1 ) * nbins_aerosol + 1 
2823          endi = index_no * nbins_aerosol
2824          ic = 1
2825          DO ss = str, endi
2826             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhohno3
2827             ic = ic+1
2828          ENDDO
2829          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2830       ENDIF
2831!
2832!--    Ammonium (NH(4+)) or ammonia NH3
2833       IF ( index_nh > 0 )  THEN
2834          vc = 7
2835          str = ( index_nh-1 ) * nbins_aerosol + 1
2836          endi = index_nh * nbins_aerosol
2837          ic = 1
2838          DO ss = str, endi
2839             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhonh3
2840             ic = ic+1
2841          ENDDO
2842          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2843       ENDIF
2844!
2845!--    Water (always used)
2846       nc_h2o = get_index( prtcl,'H2O' )
2847       vc = 8
2848       str = ( nc_h2o-1 ) * nbins_aerosol + 1
2849       endi = nc_h2o * nbins_aerosol
2850       ic = 1
2851       IF ( advect_particle_water )  THEN
2852          DO ss = str, endi
2853             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2o
2854             ic = ic+1
2855          ENDDO
2856       ELSE
2857         lo_aero(1:nbins_aerosol)%volc(vc) = mclim
2858       ENDIF
2859       aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2860!
2861!--    Number concentrations (numc) and particle sizes
2862!--    (dwet = wet diameter, core = dry volume)
2863       DO  ib = 1, nbins_aerosol
2864          lo_aero(ib)%numc = aerosol_number(ib)%conc(k,j,i)
2865          aero_old(ib)%numc = lo_aero(ib)%numc
2866          IF ( lo_aero(ib)%numc > nclim )  THEN
2867             lo_aero(ib)%dwet = ( SUM( lo_aero(ib)%volc(:) ) / lo_aero(ib)%numc / api6 )**0.33333333_wp
2868             lo_aero(ib)%core = SUM( lo_aero(ib)%volc(1:7) ) / lo_aero(ib)%numc
2869          ELSE
2870             lo_aero(ib)%dwet = lo_aero(ib)%dmid
2871             lo_aero(ib)%core = api6 * ( lo_aero(ib)%dwet )**3
2872          ENDIF
2873       ENDDO
2874!
2875!--    Calculate the ambient sizes of particles by equilibrating soluble fraction of particles with
2876!--    water using the ZSR method.
2877       in_rh = in_cw(k) / in_cs(k)
2878       IF ( prunmode==1  .OR.  .NOT. advect_particle_water )  THEN
2879          CALL equilibration( in_rh, in_t(k), lo_aero, .TRUE. )
2880       ENDIF
2881!
2882!--    Gaseous tracer concentrations in #/m3
2883       IF ( salsa_gases_from_chem )  THEN
2884!
2885!--       Convert concentrations in ppm to #/m3
2886          zgso4  = chem_species(gas_index_chem(1))%conc(k,j,i) * ppm_to_nconc(k)
2887          zghno3 = chem_species(gas_index_chem(2))%conc(k,j,i) * ppm_to_nconc(k)
2888          zgnh3  = chem_species(gas_index_chem(3))%conc(k,j,i) * ppm_to_nconc(k)
2889          zgocnv = chem_species(gas_index_chem(4))%conc(k,j,i) * ppm_to_nconc(k)
2890          zgocsv = chem_species(gas_index_chem(5))%conc(k,j,i) * ppm_to_nconc(k)
2891       ELSE
2892          zgso4  = salsa_gas(1)%conc(k,j,i)
2893          zghno3 = salsa_gas(2)%conc(k,j,i)
2894          zgnh3  = salsa_gas(3)%conc(k,j,i)
2895          zgocnv = salsa_gas(4)%conc(k,j,i)
2896          zgocsv = salsa_gas(5)%conc(k,j,i)
2897       ENDIF
2898!
2899!--    Calculate aerosol processes:
2900!--    *********************************************************************************************
2901!
2902!--    Coagulation
2903       IF ( lscoag )   THEN
2904          CALL coagulation( lo_aero, dt_salsa, in_t(k), in_p(k) )
2905       ENDIF
2906!
2907!--    Condensation
2908       IF ( lscnd )   THEN
2909          CALL condensation( lo_aero, zgso4, zgocnv, zgocsv,  zghno3, zgnh3, in_cw(k), in_cs(k),      &
2910                             in_t(k), in_p(k), dt_salsa, prtcl )
2911       ENDIF
2912!
2913!--    Deposition
2914       IF ( lsdepo )  THEN
2915          CALL deposition( lo_aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:),   &
2916                           vd(k,:) )
2917       ENDIF
2918!
2919!--    Size distribution bin update
2920       IF ( lsdistupdate )   THEN
2921          CALL distr_update( lo_aero )
2922       ENDIF
2923!--    *********************************************************************************************
2924
2925       IF ( lsdepo ) sedim_vd(k,j,i,:) = vd(k,:)
2926!
2927!--    Calculate changes in concentrations
2928       DO ib = 1, nbins_aerosol
2929          aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc -   &
2930                                           aero_old(ib)%numc ) * flag
2931       ENDDO
2932
2933       IF ( index_so4 > 0 )  THEN
2934          vc = 1
2935          str = ( index_so4-1 ) * nbins_aerosol + 1
2936          endi = index_so4 * nbins_aerosol
2937          ic = 1
2938          DO ss = str, endi
2939             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
2940                                            aero_old(ic)%volc(vc) ) * arhoh2so4 * flag
2941             ic = ic+1
2942          ENDDO
2943       ENDIF
2944
2945       IF ( index_oc > 0 )  THEN
2946          vc = 2
2947          str = ( index_oc-1 ) * nbins_aerosol + 1
2948          endi = index_oc * nbins_aerosol
2949          ic = 1
2950          DO ss = str, endi
2951             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
2952                                            aero_old(ic)%volc(vc) ) * arhooc * flag
2953             ic = ic+1
2954          ENDDO
2955       ENDIF
2956
2957       IF ( index_bc > 0 )  THEN
2958          vc = 3
2959          str = ( index_bc-1 ) * nbins_aerosol + 1 + end_subrange_1a
2960          endi = index_bc * nbins_aerosol
2961          ic = 1 + end_subrange_1a
2962          DO ss = str, endi
2963             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
2964                                            aero_old(ic)%volc(vc) ) * arhobc * flag
2965             ic = ic+1
2966          ENDDO
2967       ENDIF
2968
2969       IF ( index_du > 0 )  THEN
2970          vc = 4
2971          str = ( index_du-1 ) * nbins_aerosol + 1 + end_subrange_1a
2972          endi = index_du * nbins_aerosol
2973          ic = 1 + end_subrange_1a
2974          DO ss = str, endi
2975             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
2976                                            aero_old(ic)%volc(vc) ) * arhodu * flag
2977             ic = ic+1
2978          ENDDO
2979       ENDIF
2980
2981       IF ( index_ss > 0 )  THEN
2982          vc = 5
2983          str = ( index_ss-1 ) * nbins_aerosol + 1 + end_subrange_1a
2984          endi = index_ss * nbins_aerosol
2985          ic = 1 + end_subrange_1a
2986          DO ss = str, endi
2987             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
2988                                            aero_old(ic)%volc(vc) ) * arhoss * flag
2989             ic = ic+1
2990          ENDDO
2991       ENDIF
2992
2993       IF ( index_no > 0 )  THEN
2994          vc = 6
2995          str = ( index_no-1 ) * nbins_aerosol + 1
2996          endi = index_no * nbins_aerosol
2997          ic = 1
2998          DO ss = str, endi
2999             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3000                                            aero_old(ic)%volc(vc) ) * arhohno3 * flag
3001             ic = ic+1
3002          ENDDO
3003       ENDIF
3004
3005       IF ( index_nh > 0 )  THEN
3006          vc = 7
3007          str = ( index_nh-1 ) * nbins_aerosol + 1
3008          endi = index_nh * nbins_aerosol
3009          ic = 1
3010          DO ss = str, endi
3011             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3012                                            aero_old(ic)%volc(vc) ) * arhonh3 * flag
3013             ic = ic+1
3014          ENDDO
3015       ENDIF
3016
3017       IF ( advect_particle_water )  THEN
3018          nc_h2o = get_index( prtcl,'H2O' )
3019          vc = 8
3020          str = ( nc_h2o-1 ) * nbins_aerosol + 1
3021          endi = nc_h2o * nbins_aerosol
3022          ic = 1
3023          DO ss = str, endi
3024             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3025                                            aero_old(ic)%volc(vc) ) * arhoh2o * flag
3026             ic = ic+1
3027          ENDDO
3028       ENDIF
3029       IF ( prunmode == 1 )  THEN
3030          nc_h2o = get_index( prtcl,'H2O' )
3031          vc = 8
3032          str = ( nc_h2o-1 ) * nbins_aerosol + 1
3033          endi = nc_h2o * nbins_aerosol
3034          ic = 1
3035          DO ss = str, endi
3036             aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
3037                                             aero_old(ic)%volc(vc) ) * arhoh2o )
3038             IF ( k == nzb+1 )  THEN
3039                aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
3040             ELSEIF ( k == nzt  )  THEN
3041                aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
3042                aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
3043             ENDIF
3044             ic = ic+1
3045          ENDDO
3046       ENDIF
3047!
3048!--    Condensation of precursor gases
3049       IF ( lscndgas )  THEN
3050          IF ( salsa_gases_from_chem )  THEN
3051!
3052!--          SO4 (or H2SO4)
3053             ig = gas_index_chem(1)
3054             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgso4 /               &
3055                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3056!
3057!--          HNO3
3058             ig = gas_index_chem(2)
3059             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zghno3 /              &
3060                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3061!
3062!--          NH3
3063             ig = gas_index_chem(3)
3064             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgnh3 /               &
3065                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3066!
3067!--          non-volatile OC
3068             ig = gas_index_chem(4)
3069             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgocnv /              &
3070                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3071!
3072!--          semi-volatile OC
3073             ig = gas_index_chem(5)
3074             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgocsv /              &
3075                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3076
3077          ELSE
3078!
3079!--          SO4 (or H2SO4)
3080             salsa_gas(1)%conc(k,j,i) = salsa_gas(1)%conc(k,j,i) + ( zgso4 -                       &
3081                                        salsa_gas(1)%conc(k,j,i) ) * flag
3082!
3083!--          HNO3
3084             salsa_gas(2)%conc(k,j,i) = salsa_gas(2)%conc(k,j,i) + ( zghno3 -                      &
3085                                        salsa_gas(2)%conc(k,j,i) ) * flag
3086!
3087!--          NH3
3088             salsa_gas(3)%conc(k,j,i) = salsa_gas(3)%conc(k,j,i) + ( zgnh3 -                       &
3089                                        salsa_gas(3)%conc(k,j,i) ) * flag
3090!
3091!--          non-volatile OC
3092             salsa_gas(4)%conc(k,j,i) = salsa_gas(4)%conc(k,j,i) + ( zgocnv -                      &
3093                                        salsa_gas(4)%conc(k,j,i) ) * flag
3094!
3095!--          semi-volatile OC
3096             salsa_gas(5)%conc(k,j,i) = salsa_gas(5)%conc(k,j,i) + ( zgocsv -                      &
3097                                        salsa_gas(5)%conc(k,j,i) ) * flag
3098          ENDIF
3099       ENDIF
3100!
3101!--    Tendency of water vapour mixing ratio is obtained from the change in RH during SALSA run.
3102!--    This releases heat and changes pt. Assumes no temperature change during SALSA run.
3103!--    q = r / (1+r), Euler method for integration
3104!
3105       IF ( feedback_to_palm )  THEN
3106          q_p(k,j,i) = q_p(k,j,i) + 1.0_wp / ( in_cw(k) * in_adn(k) + 1.0_wp )**2 *                &
3107                       ( in_cw(k) - cw_old ) * in_adn(k) * flag
3108          pt_p(k,j,i) = pt_p(k,j,i) + alv / c_p * ( in_cw(k) - cw_old ) * in_adn(k) / ( in_cw(k) / &
3109                        in_adn(k) + 1.0_wp )**2 * pt_p(k,j,i) / in_t(k) * flag
3110       ENDIF
3111
3112    ENDDO   ! k
3113
3114!
3115!-- Set surfaces and wall fluxes due to deposition
3116    IF ( lsdepo  .AND.  lsdepo_surf  .AND.  prunmode == 3 )  THEN
3117       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
3118          CALL depo_surf( i, j, surf_def_h(0), vd, schmidt_num, kvis, in_u, .TRUE. )
3119          DO  l = 0, 3
3120             CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE. )
3121          ENDDO
3122       ELSE
3123          CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE., usm_to_depo_h )
3124          DO  l = 0, 3
3125             CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
3126                             usm_to_depo_v(l) )
3127          ENDDO
3128          CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE., lsm_to_depo_h )
3129          DO  l = 0, 3
3130             CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
3131                             lsm_to_depo_v(l) )
3132          ENDDO
3133       ENDIF
3134    ENDIF
3135
3136    IF ( prunmode < 3 )  THEN
3137       !$OMP MASTER
3138       aero = lo_aero
3139       !$OMP END MASTER
3140    END IF
3141
3142 END SUBROUTINE salsa_driver
3143
3144!------------------------------------------------------------------------------!
3145! Description:
3146! ------------
3147!> Set logical switches according to the host model state and user-specified
3148!> NAMELIST options.
3149!> Juha Tonttila, FMI, 2014
3150!> Only aerosol processes included, Mona Kurppa, UHel, 2017
3151!------------------------------------------------------------------------------!
3152 SUBROUTINE set_salsa_runtime( prunmode )
3153
3154    IMPLICIT NONE
3155
3156    INTEGER(iwp), INTENT(in) ::  prunmode
3157
3158    SELECT CASE(prunmode)
3159
3160       CASE(1) !< Initialization
3161          lscoag       = .FALSE.
3162          lscnd        = .FALSE.
3163          lscndgas     = .FALSE.
3164          lscndh2oae   = .FALSE.
3165          lsdepo       = .FALSE.
3166          lsdepo_pcm   = .FALSE.
3167          lsdepo_surf  = .FALSE.
3168          lsdistupdate = .TRUE.
3169          lspartition  = .FALSE.
3170
3171       CASE(2)  !< Spinup period
3172          lscoag      = ( .FALSE. .AND. nlcoag   )
3173          lscnd       = ( .TRUE.  .AND. nlcnd    )
3174          lscndgas    = ( .TRUE.  .AND. nlcndgas )
3175          lscndh2oae  = ( .TRUE.  .AND. nlcndh2oae )
3176
3177       CASE(3)  !< Run
3178          lscoag       = nlcoag
3179          lscnd        = nlcnd
3180          lscndgas     = nlcndgas
3181          lscndh2oae   = nlcndh2oae
3182          lsdepo       = nldepo
3183          lsdepo_pcm   = nldepo_pcm
3184          lsdepo_surf  = nldepo_surf
3185          lsdistupdate = nldistupdate
3186    END SELECT
3187
3188
3189 END SUBROUTINE set_salsa_runtime
3190 
3191!------------------------------------------------------------------------------!
3192! Description:
3193! ------------
3194!> Calculates the absolute temperature (using hydrostatic pressure), saturation
3195!> vapour pressure and mixing ratio over water, relative humidity and air
3196!> density needed in the SALSA model.
3197!> NOTE, no saturation adjustment takes place -> the resulting water vapour
3198!> mixing ratio can be supersaturated, allowing the microphysical calculations
3199!> in SALSA.
3200!
3201!> Juha Tonttila, FMI, 2014 (original SALSAthrm)
3202!> Mona Kurppa, UHel, 2017 (adjustment for PALM and only aerosol processes)
3203!------------------------------------------------------------------------------!
3204 SUBROUTINE salsa_thrm_ij( i, j, p_ij, temp_ij, cw_ij, cs_ij, adn_ij )
3205
3206    USE arrays_3d,                                                                                 &
3207        ONLY: pt, q, zu
3208
3209    USE basic_constants_and_equations_mod,                                                         &
3210        ONLY:  barometric_formula, exner_function, ideal_gas_law_rho, magnus
3211
3212    USE control_parameters,                                                                        &
3213        ONLY: pt_surface, surface_pressure
3214
3215    IMPLICIT NONE
3216
3217    INTEGER(iwp), INTENT(in) ::  i  !<
3218    INTEGER(iwp), INTENT(in) ::  j  !<
3219
3220    REAL(wp) ::  t_surface  !< absolute surface temperature (K)
3221
3222    REAL(wp), DIMENSION(nzb:nzt+1) ::  e_s  !< saturation vapour pressure over water (Pa)
3223
3224    REAL(wp), DIMENSION(:), INTENT(inout) ::  adn_ij   !< air density (kg/m3)
3225    REAL(wp), DIMENSION(:), INTENT(inout) ::  p_ij     !< air pressure (Pa)
3226    REAL(wp), DIMENSION(:), INTENT(inout) ::  temp_ij  !< air temperature (K)
3227
3228    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cw_ij  !< water vapour concentration (kg/m3)
3229    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cs_ij  !< saturation water vap. conc.(kg/m3)
3230!
3231!-- Pressure p_ijk (Pa) = hydrostatic pressure
3232    t_surface = pt_surface * exner_function( surface_pressure * 100.0_wp )
3233    p_ij(:) = barometric_formula( zu, t_surface, surface_pressure * 100.0_wp )
3234!
3235!-- Absolute ambient temperature (K)
3236    temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) )
3237!
3238!-- Air density
3239    adn_ij(:) = ideal_gas_law_rho( p_ij(:), temp_ij(:) )
3240!
3241!-- Water vapour concentration r_v (kg/m3)
3242    IF ( PRESENT( cw_ij ) )  THEN
3243       cw_ij(:) = ( q(:,j,i) / ( 1.0_wp - q(:,j,i) ) ) * adn_ij(:)
3244    ENDIF
3245!
3246!-- Saturation mixing ratio r_s (kg/kg) from vapour pressure at temp (Pa)
3247    IF ( PRESENT( cs_ij ) )  THEN
3248       e_s(:) = 611.0_wp * EXP( alv_d_rv * ( 3.6609E-3_wp - 1.0_wp /           &
3249                temp_ij(:) ) )! magnus( temp_ij(:) )
3250       cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:) - e_s(:) ) ) * adn_ij(:)
3251    ENDIF
3252
3253 END SUBROUTINE salsa_thrm_ij
3254
3255!------------------------------------------------------------------------------!
3256! Description:
3257! ------------
3258!> Calculates ambient sizes of particles by equilibrating soluble fraction of
3259!> particles with water using the ZSR method (Stokes and Robinson, 1966).
3260!> Method:
3261!> Following chemical components are assumed water-soluble
3262!> - (ammonium) sulphate (100%)
3263!> - sea salt (100 %)
3264!> - organic carbon (epsoc * 100%)
3265!> Exact thermodynamic considerations neglected.
3266!> - If particles contain no sea salt, calculation according to sulphate
3267!>   properties
3268!> - If contain sea salt but no sulphate, calculation according to sea salt
3269!>   properties
3270!> - If contain both sulphate and sea salt -> the molar fraction of these
3271!>   compounds determines which one of them is used as the basis of calculation.
3272!> If sulphate and sea salt coexist in a particle, it is assumed that the Cl is
3273!> replaced by sulphate; thus only either sulphate + organics or sea salt +
3274!> organics is included in the calculation of soluble fraction.
3275!> Molality parameterizations taken from Table 1 of Tang: Thermodynamic and
3276!> optical properties of mixed-salt aerosols of atmospheric importance,
3277!> J. Geophys. Res., 102 (D2), 1883-1893 (1997)
3278!
3279!> Coded by:
3280!> Hannele Korhonen (FMI) 2005
3281!> Harri Kokkola (FMI) 2006
3282!> Matti Niskanen(FMI) 2012
3283!> Anton Laakso  (FMI) 2013
3284!> Modified for the new aerosol datatype, Juha Tonttila (FMI) 2014
3285!
3286!> fxm: should sea salt form a solid particle when prh is very low (even though
3287!> it could be mixed with e.g. sulphate)?
3288!> fxm: crashes if no sulphate or sea salt
3289!> fxm: do we really need to consider Kelvin effect for subrange 2
3290!------------------------------------------------------------------------------!
3291 SUBROUTINE equilibration( prh, ptemp, paero, init )
3292
3293    IMPLICIT NONE
3294
3295    INTEGER(iwp) :: ib      !< loop index
3296    INTEGER(iwp) :: counti  !< loop index
3297
3298    LOGICAL, INTENT(in) ::  init   !< TRUE: Initialization, FALSE: Normal runtime: update water
3299                                   !< content only for 1a
3300
3301    REAL(wp) ::  zaw      !< water activity [0-1]
3302    REAL(wp) ::  zcore    !< Volume of dry particle
3303    REAL(wp) ::  zdold    !< Old diameter
3304    REAL(wp) ::  zdwet    !< Wet diameter or mean droplet diameter
3305    REAL(wp) ::  zke      !< Kelvin term in the Köhler equation
3306    REAL(wp) ::  zlwc     !< liquid water content [kg/m3-air]
3307    REAL(wp) ::  zrh      !< Relative humidity
3308
3309    REAL(wp), DIMENSION(maxspec) ::  zbinmol  !< binary molality of each components (mol/kg)
3310    REAL(wp), DIMENSION(maxspec) ::  zvpart   !< volume of chem. compounds in one particle
3311
3312    REAL(wp), INTENT(in) ::  prh    !< relative humidity [0-1]
3313    REAL(wp), INTENT(in) ::  ptemp  !< temperature (K)
3314
3315    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< aerosol properties
3316
3317    zaw       = 0.0_wp
3318    zlwc      = 0.0_wp
3319!
3320!-- Relative humidity:
3321    zrh = prh
3322    zrh = MAX( zrh, 0.05_wp )
3323    zrh = MIN( zrh, 0.98_wp)
3324!
3325!-- 1) Regime 1: sulphate and partly water-soluble OC. Done for every CALL
3326    DO  ib = start_subrange_1a, end_subrange_1a   ! size bin
3327
3328       zbinmol = 0.0_wp
3329       zdold   = 1.0_wp
3330       zke     = 1.02_wp
3331
3332       IF ( paero(ib)%numc > nclim )  THEN
3333!
3334!--       Volume in one particle
3335          zvpart = 0.0_wp
3336          zvpart(1:2) = paero(ib)%volc(1:2) / paero(ib)%numc
3337          zvpart(6:7) = paero(ib)%volc(6:7) / paero(ib)%numc
3338!
3339!--       Total volume and wet diameter of one dry particle
3340          zcore = SUM( zvpart(1:2) )
3341          zdwet = paero(ib)%dwet
3342
3343          counti = 0
3344          DO  WHILE ( ABS( zdwet / zdold - 1.0_wp ) > 1.0E-2_wp )
3345
3346             zdold = MAX( zdwet, 1.0E-20_wp )
3347             zaw = MAX( 1.0E-3_wp, zrh / zke ) ! To avoid underflow
3348!
3349!--          Binary molalities (mol/kg):
3350!--          Sulphate
3351             zbinmol(1) = 1.1065495E+2_wp - 3.6759197E+2_wp * zaw + 5.0462934E+2_wp * zaw**2 -     &
3352                          3.1543839E+2_wp * zaw**3 + 6.770824E+1_wp  * zaw**4
3353!--          Organic carbon
3354             zbinmol(2) = 1.0_wp / ( zaw * amh2o ) - 1.0_wp / amh2o
3355!--          Nitric acid
3356             zbinmol(6) = 2.306844303E+1_wp - 3.563608869E+1_wp * zaw - 6.210577919E+1_wp * zaw**2 &
3357                          + 5.510176187E+2_wp * zaw**3 - 1.460055286E+3_wp * zaw**4                &
3358                          + 1.894467542E+3_wp * zaw**5 - 1.220611402E+3_wp * zaw**6                &
3359                          + 3.098597737E+2_wp * zaw**7
3360!
3361!--          Calculate the liquid water content (kg/m3-air) using ZSR (see e.g. Eq. 10.98 in
3362!--          Seinfeld and Pandis (2006))
3363             zlwc = ( paero(ib)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / zbinmol(1) +                 &
3364                    epsoc * paero(ib)%volc(2) * ( arhooc / amoc ) / zbinmol(2) +                   &
3365                    ( paero(ib)%volc(6) * ( arhohno3/amhno3 ) ) / zbinmol(6)
3366!
3367!--          Particle wet diameter (m)
3368             zdwet = ( zlwc / paero(ib)%numc / arhoh2o / api6 + ( SUM( zvpart(6:7) ) / api6 ) +    &
3369                       zcore / api6 )**0.33333333_wp
3370!
3371!--          Kelvin effect (Eq. 10.85 in in Seinfeld and Pandis (2006)). Avoid
3372!--          overflow.
3373             zke = EXP( MIN( 50.0_wp, 4.0_wp * surfw0 * amvh2so4 / ( abo * ptemp *  zdwet ) ) )
3374
3375             counti = counti + 1
3376             IF ( counti > 1000 )  THEN
3377                message_string = 'Subrange 1: no convergence!'
3378                CALL message( 'salsa_mod: equilibration', 'PA0617', 1, 2, 0, 6, 0 )
3379             ENDIF
3380          ENDDO
3381!
3382!--       Instead of lwc, use the volume concentration of water from now on
3383!--       (easy to convert...)
3384          paero(ib)%volc(8) = zlwc / arhoh2o
3385!
3386!--       If this is initialization, update the core and wet diameter
3387          IF ( init )  THEN
3388             paero(ib)%dwet = zdwet
3389             paero(ib)%core = zcore
3390          ENDIF
3391
3392       ELSE
3393!--       If initialization
3394!--       1.2) empty bins given bin average values
3395          IF ( init )  THEN
3396             paero(ib)%dwet = paero(ib)%dmid
3397             paero(ib)%core = api6 * paero(ib)%dmid**3
3398          ENDIF
3399
3400       ENDIF
3401
3402    ENDDO  ! ib
3403!
3404!-- 2) Regime 2a: sulphate, OC, BC and sea salt
3405!--    This is done only for initialization call, otherwise the water contents
3406!--    are computed via condensation
3407    IF ( init )  THEN
3408       DO  ib = start_subrange_2a, end_subrange_2b
3409!
3410!--       Initialize
3411          zke     = 1.02_wp
3412          zbinmol = 0.0_wp
3413          zdold   = 1.0_wp
3414!
3415!--       1) Particle properties calculated for non-empty bins
3416          IF ( paero(ib)%numc > nclim )  THEN
3417!
3418!--          Volume in one particle [fxm]
3419             zvpart = 0.0_wp
3420             zvpart(1:7) = paero(ib)%volc(1:7) / paero(ib)%numc
3421!
3422!--          Total volume and wet diameter of one dry particle [fxm]
3423             zcore = SUM( zvpart(1:5) )
3424             zdwet = paero(ib)%dwet
3425
3426             counti = 0
3427             DO  WHILE ( ABS( zdwet / zdold - 1.0_wp ) > 1.0E-12_wp )
3428
3429                zdold = MAX( zdwet, 1.0E-20_wp )
3430                zaw = zrh / zke
3431!
3432!--             Binary molalities (mol/kg):
3433!--             Sulphate
3434                zbinmol(1) = 1.1065495E+2_wp - 3.6759197E+2_wp * zaw + 5.0462934E+2_wp * zaw**2 -  &
3435                             3.1543839E+2_wp * zaw**3 + 6.770824E+1_wp  * zaw**4
3436!--             Organic carbon
3437                zbinmol(2) = 1.0_wp / ( zaw * amh2o ) - 1.0_wp / amh2o
3438!--             Nitric acid
3439                zbinmol(6) = 2.306844303E+1_wp          - 3.563608869E+1_wp * zaw -                &
3440                             6.210577919E+1_wp * zaw**2 + 5.510176187E+2_wp * zaw**3 -             &
3441                             1.460055286E+3_wp * zaw**4 + 1.894467542E+3_wp * zaw**5 -             &
3442                             1.220611402E+3_wp * zaw**6 + 3.098597737E+2_wp * zaw**7 
3443!--             Sea salt (natrium chloride)
3444                zbinmol(5) = 5.875248E+1_wp - 1.8781997E+2_wp * zaw + 2.7211377E+2_wp * zaw**2 -   &
3445                             1.8458287E+2_wp * zaw**3 + 4.153689E+1_wp  * zaw**4
3446!
3447!--             Calculate the liquid water content (kg/m3-air)
3448                zlwc = ( paero(ib)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / zbinmol(1) +              &
3449                       epsoc * ( paero(ib)%volc(2) * ( arhooc / amoc ) ) / zbinmol(2) +            &
3450                       ( paero(ib)%volc(6) * ( arhohno3 / amhno3 ) ) / zbinmol(6) +                &
3451                       ( paero(ib)%volc(5) * ( arhoss / amss ) ) / zbinmol(5)
3452
3453!--             Particle wet radius (m)
3454                zdwet = ( zlwc / paero(ib)%numc / arhoh2o / api6 + ( SUM( zvpart(6:7) ) / api6 )  + &
3455                           zcore / api6 )**0.33333333_wp
3456!
3457!--             Kelvin effect (Eq. 10.85 in Seinfeld and Pandis (2006))
3458                zke = EXP( MIN( 50.0_wp, 4.0_wp * surfw0 * amvh2so4 / ( abo * zdwet * ptemp ) ) )
3459
3460                counti = counti + 1
3461                IF ( counti > 1000 )  THEN
3462                   message_string = 'Subrange 2: no convergence!'
3463                CALL message( 'salsa_mod: equilibration', 'PA0618', 1, 2, 0, 6, 0 )
3464                ENDIF
3465             ENDDO
3466!
3467!--          Liquid water content; instead of LWC use the volume concentration
3468             paero(ib)%volc(8) = zlwc / arhoh2o
3469             paero(ib)%dwet    = zdwet
3470             paero(ib)%core    = zcore
3471
3472          ELSE
3473!--          2.2) empty bins given bin average values
3474             paero(ib)%dwet = paero(ib)%dmid
3475             paero(ib)%core = api6 * paero(ib)%dmid**3
3476          ENDIF
3477
3478       ENDDO   ! ib
3479    ENDIF
3480
3481 END SUBROUTINE equilibration
3482
3483!------------------------------------------------------------------------------!
3484!> Description:
3485!> ------------
3486!> Calculation of the settling velocity vc (m/s) per aerosol size bin and
3487!> deposition on plant canopy (lsdepo_pcm).
3488!
3489!> Deposition is based on either the scheme presented in:
3490!> Zhang et al. (2001), Atmos. Environ. 35, 549-560 (includes collection due to
3491!> Brownian diffusion, impaction, interception and sedimentation; hereafter ZO1)
3492!> OR
3493!> Petroff & Zhang (2010), Geosci. Model Dev. 3, 753-769 (includes also
3494!> collection due to turbulent impaction, hereafter P10)
3495!
3496!> Equation numbers refer to equation in Jacobson (2005): Fundamentals of
3497!> Atmospheric Modeling, 2nd Edition.
3498!
3499!> Subroutine follows closely sedim_SALSA in UCLALES-SALSA written by Juha
3500!> Tonttila (KIT/FMI) and Zubair Maalick (UEF).
3501!> Rewritten to PALM by Mona Kurppa (UH), 2017.
3502!
3503!> Call for grid point i,j,k
3504!------------------------------------------------------------------------------!
3505
3506 SUBROUTINE deposition( paero, tk, adn, mag_u, lad, kvis, schmidt_num, vc )
3507
3508    USE plant_canopy_model_mod,                                                &
3509        ONLY: cdc
3510
3511    IMPLICIT NONE
3512
3513    INTEGER(iwp) ::  ib   !< loop index
3514    INTEGER(iwp) ::  ic   !< loop index
3515
3516    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
3517    REAL(wp) ::  avis              !< molecular viscocity of air (kg/(m*s))
3518    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
3519    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
3520    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
3521    REAL(wp) ::  c_interception    !< coefficient for interception
3522    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
3523    REAL(wp) ::  depo              !< deposition velocity (m/s)
3524    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
3525    REAL(wp) ::  lambda            !< molecular mean free path (m)
3526    REAL(wp) ::  mdiff             !< particle diffusivity coefficient
3527    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
3528                                   !< Table 3 in Z01
3529    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
3530    REAL(wp) ::  pdn               !< particle density (kg/m3)
3531    REAL(wp) ::  ustar             !< friction velocity (m/s)
3532    REAL(wp) ::  va                !< thermal speed of an air molecule (m/s)
3533
3534    REAL(wp), INTENT(in) ::  adn    !< air density (kg/m3)
3535    REAL(wp), INTENT(in) ::  lad    !< leaf area density (m2/m3)
3536    REAL(wp), INTENT(in) ::  mag_u  !< wind velocity (m/s)
3537    REAL(wp), INTENT(in) ::  tk     !< abs.temperature (K)
3538
3539    REAL(wp), INTENT(inout) ::  kvis   !< kinematic viscosity of air (m2/s)
3540
3541    REAL(wp), DIMENSION(nbins_aerosol) ::  beta   !< Cunningham slip-flow correction factor
3542    REAL(wp), DIMENSION(nbins_aerosol) ::  Kn     !< Knudsen number
3543    REAL(wp), DIMENSION(nbins_aerosol) ::  zdwet  !< wet diameter (m)
3544
3545    REAL(wp), DIMENSION(:), INTENT(inout) ::  schmidt_num  !< particle Schmidt number
3546    REAL(wp), DIMENSION(:), INTENT(inout) ::  vc  !< critical fall speed i.e. settling velocity of
3547                                                  !< an aerosol particle (m/s)
3548
3549    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< aerosol properties
3550!
3551!-- Initialise
3552    depo  = 0.0_wp
3553    pdn   = 1500.0_wp    ! default value
3554    ustar = 0.0_wp
3555!
3556!-- Molecular viscosity of air (Eq. 4.54)
3557    avis = 1.8325E-5_wp * ( 416.16_wp / ( tk + 120.0_wp ) ) * ( tk / 296.16_wp )**1.5_wp
3558!
3559!-- Kinematic viscosity (Eq. 4.55)
3560    kvis =  avis / adn
3561!
3562!-- Thermal velocity of an air molecule (Eq. 15.32)
3563    va = SQRT( 8.0_wp * abo * tk / ( pi * am_airmol ) )
3564!
3565!-- Mean free path (m) (Eq. 15.24)
3566    lambda = 2.0_wp * avis / ( adn * va )
3567!
3568!-- Particle wet diameter (m)
3569    zdwet = paero(:)%dwet
3570!
3571!-- Knudsen number (Eq. 15.23)
3572    Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
3573!
3574!-- Cunningham slip-flow correction (Eq. 15.30)
3575    beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
3576!
3577!-- Critical fall speed i.e. settling velocity  (Eq. 20.4)
3578    vc = MIN( 1.0_wp, zdwet**2 * ( pdn - adn ) * g * beta / ( 18.0_wp * avis ) )
3579!
3580!-- Deposition on vegetation
3581    IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
3582!
3583!--    Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
3584       alpha   = alpha_z01(depo_pcm_type_num)
3585       gamma   = gamma_z01(depo_pcm_type_num)
3586       par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
3587!
3588!--    Deposition efficiencies from Table 1. Constants from Table 2.
3589       par_l            = l_p10(depo_pcm_type_num) * 0.01_wp
3590       c_brownian_diff  = c_b_p10(depo_pcm_type_num)
3591       c_interception   = c_in_p10(depo_pcm_type_num)
3592       c_impaction      = c_im_p10(depo_pcm_type_num)
3593       beta_im          = beta_im_p10(depo_pcm_type_num)
3594       c_turb_impaction = c_it_p10(depo_pcm_type_num)
3595
3596       DO  ib = 1, nbins_aerosol
3597
3598          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
3599
3600!--       Particle diffusivity coefficient (Eq. 15.29)
3601          mdiff = ( abo * tk * beta(ib) ) / ( 3.0_wp * pi * avis * zdwet(ib) )
3602!
3603!--       Particle Schmidt number (Eq. 15.36)
3604          schmidt_num(ib) = kvis / mdiff
3605!
3606!--       Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
3607          ustar = SQRT( cdc ) * mag_u
3608          SELECT CASE ( depo_pcm_par_num )
3609
3610             CASE ( 1 )   ! Zhang et al. (2001)
3611                CALL depo_vel_Z01( vc(ib), ustar, schmidt_num(ib), paero(ib)%dwet, alpha,  gamma,  &
3612                                   par_a, depo )
3613             CASE ( 2 )   ! Petroff & Zhang (2010)
3614                CALL depo_vel_P10( vc(ib), mag_u, ustar, kvis, schmidt_num(ib), paero(ib)%dwet,    &
3615                                   par_l, c_brownian_diff, c_interception, c_impaction, beta_im,   &
3616                                   c_turb_impaction, depo )
3617          END SELECT
3618!
3619!--       Calculate the change in concentrations
3620          paero(ib)%numc = paero(ib)%numc - depo * lad * paero(ib)%numc * dt_salsa
3621          DO  ic = 1, maxspec+1
3622             paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa
3623          ENDDO
3624       ENDDO
3625
3626    ENDIF
3627
3628 END SUBROUTINE deposition
3629
3630!------------------------------------------------------------------------------!
3631! Description:
3632! ------------
3633!> Calculate deposition velocity (m/s) based on Zhan et al. (2001, case 1).
3634!------------------------------------------------------------------------------!
3635
3636 SUBROUTINE depo_vel_Z01( vc, ustar, schmidt_num, diameter, alpha, gamma, par_a, depo )
3637
3638    IMPLICIT NONE
3639
3640    REAL(wp) ::  rs                !< overall quasi-laminar resistance for particles
3641    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
3642
3643    REAL(wp), INTENT(in) ::  alpha        !< parameter, Table 3 in Z01
3644    REAL(wp), INTENT(in) ::  gamma        !< parameter, Table 3 in Z01
3645    REAL(wp), INTENT(in) ::  par_a        !< parameter A for the characteristic diameter of
3646                                          !< collectors, Table 3 in Z01
3647    REAL(wp), INTENT(in) ::  diameter     !< particle diameter
3648    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
3649    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
3650    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
3651
3652    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
3653
3654    IF ( par_a > 0.0_wp )  THEN
3655!
3656!--    Initialise
3657       rs = 0.0_wp
3658!
3659!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
3660       stokes_num = vc * ustar / ( g * par_a )
3661!
3662!--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)
3663       rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -stokes_num**0.5_wp ) *                &
3664                 ( schmidt_num**( -gamma ) + ( stokes_num / ( alpha + stokes_num ) )**2 +          &
3665                 0.5_wp * ( diameter / par_a )**2 ) ) )
3666
3667       depo = rs + vc
3668
3669    ELSE
3670       depo = 0.0_wp
3671    ENDIF
3672
3673 END SUBROUTINE depo_vel_Z01
3674
3675!------------------------------------------------------------------------------!
3676! Description:
3677! ------------
3678!> Calculate deposition velocity (m/s) based on Petroff & Zhang (2010, case 2).
3679!------------------------------------------------------------------------------!
3680
3681 SUBROUTINE depo_vel_P10( vc, mag_u, ustar, kvis_a, schmidt_num, diameter, par_l, c_brownian_diff, &
3682                          c_interception, c_impaction, beta_im, c_turb_impaction, depo )
3683
3684    IMPLICIT NONE
3685
3686    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
3687    REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
3688    REAL(wp) ::  v_bd              !< deposition velocity due to Brownian diffusion
3689    REAL(wp) ::  v_im              !< deposition velocity due to impaction
3690    REAL(wp) ::  v_in              !< deposition velocity due to interception
3691    REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
3692
3693    REAL(wp), INTENT(in) ::  beta_im           !< parameter for turbulent impaction
3694    REAL(wp), INTENT(in) ::  c_brownian_diff   !< coefficient for Brownian diffusion
3695    REAL(wp), INTENT(in) ::  c_impaction       !< coefficient for inertial impaction
3696    REAL(wp), INTENT(in) ::  c_interception    !< coefficient for interception
3697    REAL(wp), INTENT(in) ::  c_turb_impaction  !< coefficient for turbulent impaction
3698    REAL(wp), INTENT(in) ::  kvis_a       !< kinematic viscosity of air (m2/s)
3699    REAL(wp), INTENT(in) ::  mag_u        !< wind velocity (m/s)
3700    REAL(wp), INTENT(in) ::  par_l        !< obstacle characteristic dimension in P10
3701    REAL(wp), INTENT(in) ::  diameter       !< particle diameter
3702    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
3703    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
3704    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
3705
3706    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
3707
3708    IF ( par_l > 0.0_wp )  THEN
3709!
3710!--    Initialise
3711       tau_plus = 0.0_wp
3712       v_bd     = 0.0_wp
3713       v_im     = 0.0_wp
3714       v_in     = 0.0_wp
3715       v_it     = 0.0_wp
3716!
3717!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
3718       stokes_num = vc * ustar / ( g * par_l )
3719!
3720!--    Non-dimensional relexation time of the particle on top of canopy
3721       tau_plus = vc * ustar**2 / ( kvis_a * g )
3722!
3723!--    Brownian diffusion
3724       v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                             &
3725              ( mag_u * par_l / kvis_a )**( -0.5_wp )
3726!
3727!--    Interception
3728       v_in = mag_u * c_interception * diameter / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) )
3729!
3730!--    Impaction: Petroff (2009) Eq. 18
3731       v_im = mag_u * c_impaction * ( stokes_num / ( stokes_num + beta_im ) )**2
3732!
3733!--    Turbulent impaction
3734       IF ( tau_plus < 20.0_wp )  THEN
3735          v_it = 2.5E-3_wp * c_turb_impaction * tau_plus**2
3736       ELSE
3737          v_it = c_turb_impaction
3738       ENDIF
3739
3740       depo = ( v_bd + v_in + v_im + v_it + vc )
3741
3742    ELSE
3743       depo = 0.0_wp
3744    ENDIF
3745
3746 END SUBROUTINE depo_vel_P10
3747
3748!------------------------------------------------------------------------------!
3749! Description:
3750! ------------
3751!> Calculate the dry deposition on horizontal and vertical surfaces. Implement
3752!> as a surface flux.
3753!> @todo aerodynamic resistance ignored for now (not important for
3754!        high-resolution simulations)
3755!------------------------------------------------------------------------------!
3756 SUBROUTINE depo_surf( i, j, surf, vc, schmidt_num, kvis, mag_u, norm, match_array )
3757
3758    USE arrays_3d,                                                                                 &
3759        ONLY: rho_air_zw
3760
3761    USE surface_mod,                                                                               &
3762        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
3763
3764    IMPLICIT NONE
3765
3766    INTEGER(iwp) ::  ib      !< loop index
3767    INTEGER(iwp) ::  ic      !< loop index
3768    INTEGER(iwp) ::  icc     !< additional loop index
3769    INTEGER(iwp) ::  k       !< loop index
3770    INTEGER(iwp) ::  m       !< loop index
3771    INTEGER(iwp) ::  surf_e  !< End index of surface elements at (j,i)-gridpoint
3772    INTEGER(iwp) ::  surf_s  !< Start index of surface elements at (j,i)-gridpoint
3773
3774    INTEGER(iwp), INTENT(in) ::  i  !< loop index
3775    INTEGER(iwp), INTENT(in) ::  j  !< loop index
3776
3777    LOGICAL, INTENT(in) ::  norm   !< to normalise or not
3778
3779    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
3780    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
3781    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
3782    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
3783    REAL(wp) ::  c_interception    !< coefficient for interception
3784    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
3785    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
3786    REAL(wp) ::  norm_fac          !< normalisation factor (usually air density)
3787    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
3788                                   !< Table 3 in Z01
3789    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
3790    REAL(wp) ::  rs                !< the overall quasi-laminar resistance for particles
3791    REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
3792    REAL(wp) ::  v_bd              !< deposition velocity due to Brownian diffusion
3793    REAL(wp) ::  v_im              !< deposition velocity due to impaction
3794    REAL(wp) ::  v_in              !< deposition velocity due to interception
3795    REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
3796
3797    REAL(wp), DIMENSION(nbins_aerosol) ::  depo      !< deposition efficiency
3798    REAL(wp), DIMENSION(nbins_aerosol) ::  depo_sum  !< sum of deposition efficiencies
3799
3800    REAL(wp), DIMENSION(:), INTENT(in) ::  kvis   !< kinematic viscosity of air (m2/s)
3801    REAL(wp), DIMENSION(:), INTENT(in) ::  mag_u  !< wind velocity (m/s)
3802
3803    REAL(wp), DIMENSION(:,:), INTENT(in) ::  schmidt_num   !< particle Schmidt number
3804    REAL(wp), DIMENSION(:,:), INTENT(in) ::  vc            !< terminal velocity (m/s)
3805
3806    TYPE(match_surface), INTENT(in), OPTIONAL ::  match_array  !< match the deposition module and
3807                                                               !< LSM/USM surfaces
3808    TYPE(surf_type), INTENT(inout) :: surf                     !< respective surface type
3809!
3810!-- Initialise
3811    depo     = 0.0_wp
3812    depo_sum = 0.0_wp
3813    rs       = 0.0_wp
3814    surf_s   = surf%start_index(j,i)
3815    surf_e   = surf%end_index(j,i)
3816    tau_plus = 0.0_wp
3817    v_bd     = 0.0_wp
3818    v_im     = 0.0_wp
3819    v_in     = 0.0_wp
3820    v_it     = 0.0_wp
3821!
3822!-- Model parameters for the land use category. If LSM or USM is applied, import
3823!-- characteristics. Otherwise, apply surface type "urban".
3824    alpha   = alpha_z01(luc_urban)
3825    gamma   = gamma_z01(luc_urban)
3826    par_a   = A_z01(luc_urban, season) * 1.0E-3_wp
3827
3828    par_l            = l_p10(luc_urban) * 0.01_wp
3829    c_brownian_diff  = c_b_p10(luc_urban)
3830    c_interception   = c_in_p10(luc_urban)
3831    c_impaction      = c_im_p10(luc_urban)
3832    beta_im          = beta_im_p10(luc_urban)
3833    c_turb_impaction = c_it_p10(luc_urban)
3834
3835
3836    IF ( PRESENT( match_array ) )  THEN  ! land or urban surface model
3837
3838       DO  m = surf_s, surf_e
3839
3840          k = surf%k(m)
3841          norm_fac = 1.0_wp
3842
3843          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
3844
3845          IF ( match_array%match_lupg(m) > 0 )  THEN
3846             alpha = alpha_z01( match_array%match_lupg(m) )
3847             gamma = gamma_z01( match_array%match_lupg(m) )
3848             par_a = A_z01( match_array%match_lupg(m), season ) * 1.0E-3_wp
3849
3850             beta_im          = beta_im_p10( match_array%match_lupg(m) )
3851             c_brownian_diff  = c_b_p10( match_array%match_lupg(m) )
3852             c_impaction      = c_im_p10( match_array%match_lupg(m) )
3853             c_interception   = c_in_p10( match_array%match_lupg(m) )
3854             c_turb_impaction = c_it_p10( match_array%match_lupg(m) )
3855             par_l            = l_p10( match_array%match_lupg(m) ) * 0.01_wp
3856
3857             DO  ib = 1, nbins_aerosol
3858                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
3859                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
3860
3861                SELECT CASE ( depo_surf_par_num )
3862
3863                   CASE ( 1 )
3864                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
3865                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
3866                   CASE ( 2 )
3867                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
3868                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
3869                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
3870                                         c_turb_impaction, depo(ib) )
3871                END SELECT
3872             ENDDO
3873             depo_sum = depo_sum + surf%frac(ind_pav_green,m) * depo
3874          ENDIF
3875
3876          IF ( match_array%match_luvw(m) > 0 )  THEN
3877             alpha = alpha_z01( match_array%match_luvw(m) )
3878             gamma = gamma_z01( match_array%match_luvw(m) )
3879             par_a = A_z01( match_array%match_luvw(m), season ) * 1.0E-3_wp
3880
3881             beta_im          = beta_im_p10( match_array%match_luvw(m) )
3882             c_brownian_diff  = c_b_p10( match_array%match_luvw(m) )
3883             c_impaction      = c_im_p10( match_array%match_luvw(m) )
3884             c_interception   = c_in_p10( match_array%match_luvw(m) )
3885             c_turb_impaction = c_it_p10( match_array%match_luvw(m) )
3886             par_l            = l_p10( match_array%match_luvw(m) ) * 0.01_wp
3887
3888             DO  ib = 1, nbins_aerosol
3889                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
3890                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
3891
3892                SELECT CASE ( depo_surf_par_num )
3893
3894                   CASE ( 1 )
3895                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
3896                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
3897                   CASE ( 2 )
3898                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
3899                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
3900                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
3901                                         c_turb_impaction, depo(ib) )
3902                END SELECT
3903             ENDDO
3904             depo_sum = depo_sum + surf%frac(ind_veg_wall,m) * depo
3905          ENDIF
3906
3907          IF ( match_array%match_luww(m) > 0 )  THEN
3908             alpha = alpha_z01( match_array%match_luww(m) )
3909             gamma = gamma_z01( match_array%match_luww(m) )
3910             par_a = A_z01( match_array%match_luww(m), season ) * 1.0E-3_wp
3911
3912             beta_im          = beta_im_p10( match_array%match_luww(m) )
3913             c_brownian_diff  = c_b_p10( match_array%match_luww(m) )
3914             c_impaction      = c_im_p10( match_array%match_luww(m) )
3915             c_interception   = c_in_p10( match_array%match_luww(m) )
3916             c_turb_impaction = c_it_p10( match_array%match_luww(m) )
3917             par_l            = l_p10( match_array%match_luww(m) ) * 0.01_wp
3918
3919             DO  ib = 1, nbins_aerosol
3920                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
3921                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
3922
3923                SELECT CASE ( depo_surf_par_num )
3924
3925                   CASE ( 1 )
3926                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
3927                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
3928                   CASE ( 2 )
3929                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
3930                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
3931                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
3932                                         c_turb_impaction, depo(ib) )
3933                END SELECT
3934             ENDDO
3935             depo_sum = depo_sum + surf%frac(ind_wat_win,m) * depo
3936          ENDIF
3937
3938          DO  ib = 1, nbins_aerosol
3939             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) )  CYCLE
3940!
3941!--          Calculate changes in surface fluxes due to dry deposition
3942             IF ( include_emission )  THEN
3943                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
3944                                   depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
3945                DO  ic = 1, ncomponents_mass
3946                   icc = ( ic - 1 ) * nbins_aerosol + ib
3947                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
3948                                       depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
3949                ENDDO  ! ic
3950             ELSE
3951                surf%answs(m,ib) = -depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
3952                DO  ic = 1, ncomponents_mass
3953                   icc = ( ic - 1 ) * nbins_aerosol + ib
3954                   surf%amsws(m,icc) = -depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
3955                ENDDO  ! ic
3956             ENDIF
3957          ENDDO  ! ib
3958
3959       ENDDO
3960
3961    ELSE  ! default surfaces
3962
3963       DO  m = surf_s, surf_e
3964
3965          k = surf%k(m)
3966          norm_fac = 1.0_wp
3967
3968          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
3969
3970          DO  ib = 1, nbins_aerosol
3971             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                        &
3972                  schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
3973
3974             SELECT CASE ( depo_surf_par_num )
3975
3976                CASE ( 1 )
3977                   CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),                 &
3978                                      ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
3979                CASE ( 2 )
3980                   CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),               &
3981                                      schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,                &
3982                                      c_brownian_diff, c_interception, c_impaction, beta_im,       &
3983                                      c_turb_impaction, depo(ib) )
3984             END SELECT
3985!
3986!--          Calculate changes in surface fluxes due to dry deposition
3987             IF ( include_emission )  THEN
3988                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
3989                                   depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
3990                DO  ic = 1, ncomponents_mass
3991                   icc = ( ic - 1 ) * nbins_aerosol + ib
3992                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
3993                                       depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
3994                ENDDO  ! ic
3995             ELSE
3996                surf%answs(m,ib) = -depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
3997                DO  ic = 1, ncomponents_mass
3998                   icc = ( ic - 1 ) * nbins_aerosol + ib
3999                   surf%amsws(m,icc) = -depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
4000                ENDDO  ! ic
4001             ENDIF
4002          ENDDO  ! ib
4003       ENDDO
4004
4005    ENDIF
4006
4007 END SUBROUTINE depo_surf
4008
4009!------------------------------------------------------------------------------!
4010! Description:
4011! ------------
4012!> Calculates particle loss and change in size distribution due to (Brownian)
4013!> coagulation. Only for particles with dwet < 30 micrometres.
4014!
4015!> Method:
4016!> Semi-implicit, non-iterative method: (Jacobson, 1994)
4017!> Volume concentrations of the smaller colliding particles added to the bin of
4018!> the larger colliding particles. Start from first bin and use the updated
4019!> number and volume for calculation of following bins. NB! Our bin numbering
4020!> does not follow particle size in subrange 2.
4021!
4022!> Schematic for bin numbers in different subranges:
4023!>             1                            2
4024!>    +-------------------------------------------+
4025!>  a | 1 | 2 | 3 || 4 | 5 | 6 | 7 |  8 |  9 | 10||
4026!>  b |           ||11 |12 |13 |14 | 15 | 16 | 17||
4027!>    +-------------------------------------------+
4028!
4029!> Exact coagulation coefficients for each pressure level are scaled according
4030!> to current particle wet size (linear scaling).
4031!> Bins are organized in terms of the dry size of the condensation nucleus,
4032!> while coagulation kernell is calculated with the actual hydrometeor
4033!> size.
4034!
4035!> Called from salsa_driver
4036!> fxm: Process selection should be made smarter - now just lots of IFs inside
4037!>      loops
4038!
4039!> Coded by:
4040!> Hannele Korhonen (FMI) 2005
4041!> Harri Kokkola (FMI) 2006
4042!> Tommi Bergman (FMI) 2012
4043!> Matti Niskanen(FMI) 2012
4044!> Anton Laakso  (FMI) 2013
4045!> Juha Tonttila (FMI) 2014
4046!------------------------------------------------------------------------------!
4047 SUBROUTINE coagulation( paero, ptstep, ptemp, ppres )
4048
4049    IMPLICIT NONE
4050
4051    INTEGER(iwp) ::  index_2a !< corresponding bin in subrange 2a
4052    INTEGER(iwp) ::  index_2b !< corresponding bin in subrange 2b
4053    INTEGER(iwp) ::  ib       !< loop index
4054    INTEGER(iwp) ::  ll       !< loop index
4055    INTEGER(iwp) ::  mm       !< loop index
4056    INTEGER(iwp) ::  nn       !< loop index
4057
4058    REAL(wp) ::  pressi          !< pressure
4059    REAL(wp) ::  temppi          !< temperature
4060    REAL(wp) ::  zdpart_mm       !< diameter of particle (m)
4061    REAL(wp) ::  zdpart_nn       !< diameter of particle (m)
4062    REAL(wp) ::  zminusterm      !< coagulation loss in a bin (1/s)
4063
4064    REAL(wp), INTENT(in) ::  ppres  !< ambient pressure (Pa)
4065    REAL(wp), INTENT(in) ::  ptemp  !< ambient temperature (K)
4066    REAL(wp), INTENT(in) ::  ptstep !< time step (s)
4067
4068    REAL(wp), DIMENSION(nbins_aerosol) ::  zmpart     !< approximate mass of particles (kg)
4069    REAL(wp), DIMENSION(maxspec+1)     ::  zplusterm  !< coagulation gain in a bin (for each
4070                                                      !< chemical compound)
4071    REAL(wp), DIMENSION(nbins_aerosol,nbins_aerosol) ::  zcc  !< updated coagulation coefficients (m3/s)
4072
4073    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< Aerosol properties
4074
4075    zdpart_mm = 0.0_wp
4076    zdpart_nn = 0.0_wp
4077!
4078!-- 1) Coagulation to coarse mode calculated in a simplified way:
4079!--    CoagSink ~ Dp in continuum subrange, thus we calculate 'effective'
4080!--    number concentration of coarse particles
4081
4082!-- 2) Updating coagulation coefficients
4083!
4084!-- Aerosol mass (kg). Density of 1500 kg/m3 assumed
4085    zmpart(1:end_subrange_2b) = api6 * ( MIN( paero(1:end_subrange_2b)%dwet, 30.0E-6_wp )**3 )     &
4086                                * 1500.0_wp
4087    temppi = ptemp
4088    pressi = ppres
4089    zcc    = 0.0_wp
4090!
4091!-- Aero-aero coagulation
4092    DO  mm = 1, end_subrange_2b   ! smaller colliding particle
4093       IF ( paero(mm)%numc < ( 2.0_wp * nclim ) )  CYCLE
4094       DO  nn = mm, end_subrange_2b   ! larger colliding particle
4095          IF ( paero(nn)%numc < ( 2.0_wp * nclim ) )  CYCLE
4096
4097          zdpart_mm = MIN( paero(mm)%dwet, 30.0E-6_wp )     ! Limit to 30 um
4098          zdpart_nn = MIN( paero(nn)%dwet, 30.0E-6_wp )     ! Limit to 30 um
4099!
4100!--       Coagulation coefficient of particles (m3/s)
4101          zcc(mm,nn) = coagc( zdpart_mm, zdpart_nn, zmpart(mm), zmpart(nn), temppi, pressi )
4102          zcc(nn,mm) = zcc(mm,nn)
4103       ENDDO
4104    ENDDO
4105
4106!
4107!-- 3) New particle and volume concentrations after coagulation:
4108!--    Calculated according to Jacobson (2005) eq. 15.9
4109!
4110!-- Aerosols in subrange 1a:
4111    DO  ib = start_subrange_1a, end_subrange_1a
4112       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
4113       zminusterm   = 0.0_wp
4114       zplusterm(:) = 0.0_wp
4115!
4116!--    Particles lost by coagulation with larger aerosols
4117       DO  ll = ib+1, end_subrange_2b
4118          zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4119       ENDDO
4120!
4121!--    Coagulation gain in a bin: change in volume conc. (cm3/cm3):
4122       DO ll = start_subrange_1a, ib - 1
4123          zplusterm(1:2) = zplusterm(1:2) + zcc(ll,ib) * paero(ll)%volc(1:2)
4124          zplusterm(6:7) = zplusterm(6:7) + zcc(ll,ib) * paero(ll)%volc(6:7)
4125          zplusterm(8)   = zplusterm(8)   + zcc(ll,ib) * paero(ll)%volc(8)
4126       ENDDO
4127!
4128!--    Volume and number concentrations after coagulation update [fxm]
4129       paero(ib)%volc(1:2) = ( paero(ib)%volc(1:2) + ptstep * zplusterm(1:2) * paero(ib)%numc ) /  &
4130                            ( 1.0_wp + ptstep * zminusterm )
4131       paero(ib)%volc(6:8) = ( paero(ib)%volc(6:8) + ptstep * zplusterm(6:8) * paero(ib)%numc ) /  &
4132                            ( 1.0_wp + ptstep * zminusterm )
4133       paero(ib)%numc = paero(ib)%numc / ( 1.0_wp + ptstep * zminusterm + 0.5_wp * ptstep *        &
4134                        zcc(ib,ib) * paero(ib)%numc )
4135    ENDDO
4136!
4137!-- Aerosols in subrange 2a:
4138    DO  ib = start_subrange_2a, end_subrange_2a
4139       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
4140       zminusterm   = 0.0_wp
4141       zplusterm(:) = 0.0_wp
4142!
4143!--    Find corresponding size bin in subrange 2b
4144       index_2b = ib - start_subrange_2a + start_subrange_2b
4145!
4146!--    Particles lost by larger particles in 2a
4147       DO  ll = ib+1, end_subrange_2a
4148          zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4149       ENDDO
4150!
4151!--    Particles lost by larger particles in 2b
4152       IF ( .NOT. no_insoluble )  THEN
4153          DO  ll = index_2b+1, end_subrange_2b
4154             zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4155          ENDDO
4156       ENDIF
4157!
4158!--    Particle volume gained from smaller particles in subranges 1, 2a and 2b
4159       DO  ll = start_subrange_1a, ib-1
4160          zplusterm(1:2) = zplusterm(1:2) + zcc(ll,ib) * paero(ll)%volc(1:2)
4161          zplusterm(6:8) = zplusterm(6:8) + zcc(ll,ib) * paero(ll)%volc(6:8)
4162       ENDDO
4163!
4164!--    Particle volume gained from smaller particles in 2a
4165!--    (Note, for components not included in the previous loop!)
4166       DO  ll = start_subrange_2a, ib-1
4167          zplusterm(3:5) = zplusterm(3:5) + zcc(ll,ib)*paero(ll)%volc(3:5)
4168       ENDDO
4169!
4170!--    Particle volume gained from smaller (and equal) particles in 2b
4171       IF ( .NOT. no_insoluble )  THEN
4172          DO  ll = start_subrange_2b, index_2b
4173             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,ib) * paero(ll)%volc(1:8)
4174          ENDDO
4175       ENDIF
4176!
4177!--    Volume and number concentrations after coagulation update [fxm]
4178       paero(ib)%volc(1:8) = ( paero(ib)%volc(1:8) + ptstep * zplusterm(1:8) * paero(ib)%numc ) /  &
4179                            ( 1.0_wp + ptstep * zminusterm )
4180       paero(ib)%numc = paero(ib)%numc / ( 1.0_wp + ptstep * zminusterm + 0.5_wp * ptstep *        &
4181                        zcc(ib,ib) * paero(ib)%numc )
4182    ENDDO
4183!
4184!-- Aerosols in subrange 2b:
4185    IF ( .NOT. no_insoluble )  THEN
4186       DO  ib = start_subrange_2b, end_subrange_2b
4187          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
4188          zminusterm   = 0.0_wp
4189          zplusterm(:) = 0.0_wp
4190!
4191!--       Find corresponding size bin in subsubrange 2a
4192          index_2a = ib - start_subrange_2b + start_subrange_2a
4193!
4194!--       Particles lost to larger particles in subranges 2b
4195          DO  ll = ib + 1, end_subrange_2b
4196             zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4197          ENDDO
4198!
4199!--       Particles lost to larger and equal particles in 2a
4200          DO  ll = index_2a, end_subrange_2a
4201             zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4202          ENDDO
4203!
4204!--       Particle volume gained from smaller particles in subranges 1 & 2a
4205          DO  ll = start_subrange_1a, index_2a - 1
4206             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,ib) * paero(ll)%volc(1:8)
4207          ENDDO
4208!
4209!--       Particle volume gained from smaller particles in 2b
4210          DO  ll = start_subrange_2b, ib - 1
4211             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,ib) * paero(ll)%volc(1:8)
4212          ENDDO
4213!
4214!--       Volume and number concentrations after coagulation update [fxm]
4215          paero(ib)%volc(1:8) = ( paero(ib)%volc(1:8) + ptstep * zplusterm(1:8) * paero(ib)%numc ) &
4216                                / ( 1.0_wp + ptstep * zminusterm )
4217          paero(ib)%numc = paero(ib)%numc / ( 1.0_wp + ptstep * zminusterm + 0.5_wp * ptstep *     &
4218                           zcc(ib,ib) * paero(ib)%numc )
4219       ENDDO
4220    ENDIF
4221
4222 END SUBROUTINE coagulation
4223
4224!------------------------------------------------------------------------------!
4225! Description:
4226! ------------
4227!> Calculation of coagulation coefficients. Extended version of the function
4228!> originally found in mo_salsa_init.
4229!
4230!> J. Tonttila, FMI, 05/2014
4231!------------------------------------------------------------------------------!
4232 REAL(wp) FUNCTION coagc( diam1, diam2, mass1, mass2, temp, pres )
4233
4234    IMPLICIT NONE
4235
4236    REAL(wp) ::  fmdist  !< distance of flux matching (m)
4237    REAL(wp) ::  knud_p  !< particle Knudsen number
4238    REAL(wp) ::  mdiam   !< mean diameter of colliding particles (m)
4239    REAL(wp) ::  mfp     !< mean free path of air molecules (m)
4240    REAL(wp) ::  visc    !< viscosity of air (kg/(m s))
4241
4242    REAL(wp), INTENT(in) ::  diam1  !< diameter of colliding particle 1 (m)
4243    REAL(wp), INTENT(in) ::  diam2  !< diameter of colliding particle 2 (m)
4244    REAL(wp), INTENT(in) ::  mass1  !< mass of colliding particle 1 (kg)
4245    REAL(wp), INTENT(in) ::  mass2  !< mass of colliding particle 2 (kg)
4246    REAL(wp), INTENT(in) ::  pres   !< ambient pressure (Pa?) [fxm]
4247    REAL(wp), INTENT(in) ::  temp   !< ambient temperature (K)
4248
4249    REAL(wp), DIMENSION (2) ::  beta    !< Cunningham correction factor
4250    REAL(wp), DIMENSION (2) ::  dfpart  !< particle diffusion coefficient (m2/s)
4251    REAL(wp), DIMENSION (2) ::  diam    !< diameters of particles (m)
4252    REAL(wp), DIMENSION (2) ::  flux    !< flux in continuum and free molec. regime (m/s)
4253    REAL(wp), DIMENSION (2) ::  knud    !< particle Knudsen number
4254    REAL(wp), DIMENSION (2) ::  mpart   !< masses of particles (kg)
4255    REAL(wp), DIMENSION (2) ::  mtvel   !< particle mean thermal velocity (m/s)
4256    REAL(wp), DIMENSION (2) ::  omega   !< particle mean free path
4257    REAL(wp), DIMENSION (2) ::  tva     !< temporary variable (m)
4258!
4259!-- Initialisation
4260    coagc   = 0.0_wp
4261!
4262!-- 1) Initializing particle and ambient air variables
4263    diam  = (/ diam1, diam2 /) !< particle diameters (m)
4264    mpart = (/ mass1, mass2 /) !< particle masses (kg)
4265!
4266!-- Viscosity of air (kg/(m s))
4267    visc = ( 7.44523E-3_wp * temp ** 1.5_wp ) / ( 5093.0_wp * ( temp + 110.4_wp ) )
4268!
4269!-- Mean free path of air (m)
4270    mfp = ( 1.656E-10_wp * temp + 1.828E-8_wp ) * ( p_0 + 1325.0_wp ) / pres
4271!
4272!-- 2) Slip correction factor for small particles
4273    knud = 2.0_wp * EXP( LOG(mfp) - LOG(diam) )! Knudsen number for air (15.23)
4274!
4275!-- Cunningham correction factor (Allen and Raabe, Aerosol Sci. Tech. 4, 269)
4276    beta = 1.0_wp + knud * ( 1.142_wp + 0.558_wp * EXP( -0.999_wp / knud ) )
4277!
4278!-- 3) Particle properties
4279!-- Diffusion coefficient (m2/s) (Jacobson (2005) eq. 15.29)
4280    dfpart = beta * abo * temp / ( 3.0_wp * pi * visc * diam )
4281!
4282!-- Mean thermal velocity (m/s) (Jacobson (2005) eq. 15.32)
4283    mtvel = SQRT( ( 8.0_wp * abo * temp ) / ( pi * mpart ) )
4284!
4285!-- Particle mean free path (m) (Jacobson (2005) eq. 15.34 )
4286    omega = 8.0_wp * dfpart / ( pi * mtvel )
4287!
4288!-- Mean diameter (m)
4289    mdiam = 0.5_wp * ( diam(1) + diam(2) )
4290!
4291!-- 4) Calculation of fluxes (Brownian collision kernels) and flux matching
4292!-- following Jacobson (2005):
4293!
4294!-- Flux in continuum regime (m3/s) (eq. 15.28)
4295    flux(1) = 4.0_wp * pi * mdiam * ( dfpart(1) + dfpart(2) )
4296!
4297!-- Flux in free molec. regime (m3/s) (eq. 15.31)
4298    flux(2) = pi * SQRT( ( mtvel(1)**2 ) + ( mtvel(2)**2 ) ) * ( mdiam**2 )
4299!
4300!-- temporary variables (m) to calculate flux matching distance (m)
4301    tva(1) = ( ( mdiam + omega(1) )**3 - ( mdiam**2 + omega(1)**2 ) * SQRT( ( mdiam**2 +           &
4302               omega(1)**2 ) ) ) / ( 3.0_wp * mdiam * omega(1) ) - mdiam
4303    tva(2) = ( ( mdiam + omega(2) )**3 - ( mdiam**2 + omega(2)**2 ) * SQRT( ( mdiam**2 +           &
4304               omega(2)**2 ) ) ) / ( 3.0_wp * mdiam * omega(2) ) - mdiam
4305!
4306!-- Flux matching distance (m): the mean distance from the centre of a sphere reached by particles
4307!-- that leave sphere's surface and travel a distance of particle mean free path (eq. 15.34)
4308    fmdist = SQRT( tva(1)**2 + tva(2)**2 )
4309!
4310!-- 5) Coagulation coefficient = coalescence efficiency * collision kernel (m3/s) (eq. 15.33).
4311!--    Here assumed coalescence efficiency 1!!
4312    coagc = flux(1) / ( mdiam / ( mdiam + fmdist) + flux(1) / flux(2) )
4313!
4314!-- Corrected collision kernel (Karl et al., 2016 (ACP)): Include van der Waals and viscous forces
4315    IF ( van_der_waals_coagc )  THEN
4316       knud_p = SQRT( omega(1)**2 + omega(2)**2 ) / mdiam
4317       IF ( knud_p >= 0.1_wp  .AND.  knud_p <= 10.0_wp )  THEN
4318          coagc = coagc * ( 2.0_wp + 0.4_wp * LOG( knud_p ) )
4319       ELSE
4320          coagc = coagc * 3.0_wp
4321       ENDIF
4322    ENDIF
4323
4324 END FUNCTION coagc
4325
4326!------------------------------------------------------------------------------!
4327! Description:
4328! ------------
4329!> Calculates the change in particle volume and gas phase
4330!> concentrations due to nucleation, condensation and dissolutional growth.
4331!
4332!> Sulphuric acid and organic vapour: only condensation and no evaporation.
4333!
4334!> New gas and aerosol phase concentrations calculated according to Jacobson
4335!> (1997): Numerical techniques to solve condensational and dissolutional growth
4336!> equations when growth is coupled to reversible reactions, Aerosol Sci. Tech.,
4337!> 27, pp 491-498.
4338!
4339!> Following parameterization has been used:
4340!> Molecular diffusion coefficient of condensing vapour (m2/s)
4341!> (Reid et al. (1987): Properties of gases and liquids, McGraw-Hill, New York.)
4342!> D = {1.d-7*sqrt(1/M_air + 1/M_gas)*T^1.75} / &
4343!      {p_atm/p_stand * (d_air^(1/3) + d_gas^(1/3))^2 }
4344!> M_air = 28.965 : molar mass of air (g/mol)
4345!> d_air = 19.70  : diffusion volume of air
4346!> M_h2so4 = 98.08 : molar mass of h2so4 (g/mol)
4347!> d_h2so4 = 51.96  : diffusion volume of h2so4
4348!
4349!> Called from main aerosol model
4350!> For equations, see Jacobson, Fundamentals of Atmospheric Modeling, 2nd Edition (2005)
4351!
4352!> Coded by:
4353!> Hannele Korhonen (FMI) 2005
4354!> Harri Kokkola (FMI) 2006
4355!> Juha Tonttila (FMI) 2014
4356!> Rewritten to PALM by Mona Kurppa (UHel) 2017
4357!------------------------------------------------------------------------------!
4358 SUBROUTINE condensation( paero, pc_sa, pc_ocnv, pcocsv, pchno3, pc_nh3, pcw, pcs, ptemp, ppres,   &
4359                          ptstep, prtcl )
4360
4361    IMPLICIT NONE
4362
4363    INTEGER(iwp) ::  ss      !< start index
4364    INTEGER(iwp) ::  ee      !< end index
4365
4366    REAL(wp) ::  zcs_ocnv    !< condensation sink of nonvolatile organics (1/s)
4367    REAL(wp) ::  zcs_ocsv    !< condensation sink of semivolatile organics (1/s)
4368    REAL(wp) ::  zcs_su      !< condensation sink of sulfate (1/s)
4369    REAL(wp) ::  zcs_tot     !< total condensation sink (1/s) (gases)
4370    REAL(wp) ::  zcvap_new1  !< vapour concentration after time step (#/m3): sulphuric acid
4371    REAL(wp) ::  zcvap_new2  !< nonvolatile organics
4372    REAL(wp) ::  zcvap_new3  !< semivolatile organics
4373    REAL(wp) ::  zdfvap      !< air diffusion coefficient (m2/s)
4374    REAL(wp) ::  zdvap1      !< change in vapour concentration (#/m3): sulphuric acid
4375    REAL(wp) ::  zdvap2      !< nonvolatile organics
4376    REAL(wp) ::  zdvap3      !< semivolatile organics
4377    REAL(wp) ::  zmfp        !< mean free path of condensing vapour (m)
4378    REAL(wp) ::  zrh         !< Relative humidity [0-1]
4379    REAL(wp) ::  zvisc       !< viscosity of air (kg/(m s))
4380    REAL(wp) ::  zn_vs_c     !< ratio of nucleation of all mass transfer in the smallest bin
4381    REAL(wp) ::  zxocnv      !< ratio of organic vapour in 3nm particles
4382    REAL(wp) ::  zxsa        !< Ratio in 3nm particles: sulphuric acid
4383
4384    REAL(wp), INTENT(in) ::  ppres   !< ambient pressure (Pa)
4385    REAL(wp), INTENT(in) ::  pcs     !< Water vapour saturation concentration (kg/m3)
4386    REAL(wp), INTENT(in) ::  ptemp   !< ambient temperature (K)
4387    REAL(wp), INTENT(in) ::  ptstep  !< timestep (s)
4388
4389    REAL(wp), INTENT(inout) ::  pchno3   !< Gas concentrations (#/m3): nitric acid HNO3
4390    REAL(wp), INTENT(inout) ::  pc_nh3   !< ammonia NH3
4391    REAL(wp), INTENT(inout) ::  pc_ocnv  !< non-volatile organics
4392    REAL(wp), INTENT(inout) ::  pcocsv   !< semi-volatile organics
4393    REAL(wp), INTENT(inout) ::  pc_sa    !< sulphuric acid H2SO4
4394    REAL(wp), INTENT(inout) ::  pcw      !< Water vapor concentration (kg/m3)
4395
4396    REAL(wp), DIMENSION(nbins_aerosol)       ::  zbeta          !< transitional correction factor
4397    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate       !< collision rate (1/s)
4398    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate_ocnv  !< collision rate of OCNV (1/s)
4399    REAL(wp), DIMENSION(start_subrange_1a+1) ::  zdfpart        !< particle diffusion coef. (m2/s)
4400    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvoloc        !< change of organics volume
4401    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvolsa        !< change of sulphate volume
4402    REAL(wp), DIMENSION(2)                   ::  zj3n3          !< Formation massrate of molecules
4403                                                                !< in nucleation, (molec/m3s),
4404                                                                !< 1: H2SO4 and 2: organic vapor
4405    REAL(wp), DIMENSION(nbins_aerosol)       ::  zknud          !< particle Knudsen number
4406
4407    TYPE(component_index), INTENT(in) :: prtcl  !< Keeps track which substances are used
4408
4409    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< Aerosol properties
4410
4411    zj3n3  = 0.0_wp
4412    zrh    = pcw / pcs
4413    zxocnv = 0.0_wp
4414    zxsa   = 0.0_wp
4415!
4416!-- Nucleation
4417    IF ( nsnucl > 0 )  THEN
4418       CALL nucleation( paero, ptemp, zrh, ppres, pc_sa, pc_ocnv, pc_nh3, ptstep, zj3n3, zxsa,     &
4419                        zxocnv )
4420    ENDIF
4421!
4422!-- Condensation on pre-existing particles
4423    IF ( lscndgas )  THEN
4424!
4425!--    Initialise:
4426       zdvolsa = 0.0_wp
4427       zdvoloc = 0.0_wp
4428       zcolrate = 0.0_wp
4429!
4430!--    1) Properties of air and condensing gases:
4431!--    Viscosity of air (kg/(m s)) (Eq. 4.54 in Jabonson (2005))
4432       zvisc = ( 7.44523E-3_wp * ptemp ** 1.5_wp ) / ( 5093.0_wp * ( ptemp + 110.4_wp ) )
4433!
4434!--    Diffusion coefficient of air (m2/s)
4435       zdfvap = 5.1111E-10_wp * ptemp ** 1.75_wp * ( p_0 + 1325.0_wp ) / ppres
4436!
4437!--    Mean free path (m): same for H2SO4 and organic compounds
4438       zmfp = 3.0_wp * zdfvap * SQRT( pi * amh2so4 / ( 8.0_wp * argas * ptemp ) )
4439!
4440!--    2) Transition regime correction factor zbeta for particles (Fuchs and Sutugin (1971)):
4441!--       Size of condensing molecule considered only for nucleation mode (3 - 20 nm).
4442!
4443!--    Particle Knudsen number: condensation of gases on aerosols
4444       ss = start_subrange_1a
4445       ee = start_subrange_1a+1
4446       zknud(ss:ee) = 2.0_wp * zmfp / ( paero(ss:ee)%dwet + d_sa )
4447       ss = start_subrange_1a+2
4448       ee = end_subrange_2b
4449       zknud(ss:ee) = 2.0_wp * zmfp / paero(ss:ee)%dwet
4450!
4451!--    Transitional correction factor: aerosol + gas (the semi-empirical Fuchs- Sutugin
4452!--    interpolation function (Fuchs and Sutugin, 1971))
4453       zbeta = ( zknud + 1.0_wp ) / ( 0.377_wp * zknud + 1.0_wp + 4.0_wp / ( 3.0_wp * massacc ) *  &
4454               ( zknud + zknud ** 2 ) )
4455!
4456!--    3) Collision rate of molecules to particles
4457!--       Particle diffusion coefficient considered only for nucleation mode (3 - 20 nm)
4458!
4459!--    Particle diffusion coefficient (m2/s) (e.g. Eq. 15.29 in Jacobson (2005))
4460       zdfpart = abo * ptemp * zbeta(start_subrange_1a:start_subrange_1a+1) / ( 3.0_wp * pi * zvisc&
4461                 * paero(start_subrange_1a:start_subrange_1a+1)%dwet)
4462!
4463!--    Collision rate (mass-transfer coefficient): gases on aerosols (1/s) (Eq. 16.64 in
4464!--    Jacobson (2005))
4465       ss = start_subrange_1a
4466       ee = start_subrange_1a+1
4467       zcolrate(ss:ee) = MERGE( 2.0_wp * pi * ( paero(ss:ee)%dwet + d_sa ) * ( zdfvap + zdfpart ) *&
4468                               zbeta(ss:ee) * paero(ss:ee)%numc, 0.0_wp, paero(ss:ee)%numc > nclim )
4469       ss = start_subrange_1a+2
4470       ee = end_subrange_2b
4471       zcolrate(ss:ee) = MERGE( 2.0_wp * pi * paero(ss:ee)%dwet * zdfvap * zbeta(ss:ee) *          &
4472                                paero(ss:ee)%numc, 0.0_wp, paero(ss:ee)%numc > nclim )
4473!
4474!-- 4) Condensation sink (1/s)
4475       zcs_tot = SUM( zcolrate )   ! total sink
4476!
4477!--    5) Changes in gas-phase concentrations and particle volume
4478!
4479!--    5.1) Organic vapours
4480!
4481!--    5.1.1) Non-volatile organic compound: condenses onto all bins
4482       IF ( pc_ocnv > 1.0E+10_wp  .AND.  zcs_tot > 1.0E-30_wp  .AND. index_oc > 0 )  &
4483       THEN
4484!--       Ratio of nucleation vs. condensation rates in the smallest bin
4485          zn_vs_c = 0.0_wp
4486          IF ( zj3n3(2) > 1.0_wp )  THEN
4487             zn_vs_c = ( zj3n3(2) ) / ( zj3n3(2) + pc_ocnv * zcolrate(start_subrange_1a) )
4488          ENDIF
4489!
4490!--       Collision rate in the smallest bin, including nucleation and condensation (see
4491!--       Jacobson (2005), eq. (16.73) )
4492          zcolrate_ocnv = zcolrate
4493          zcolrate_ocnv(start_subrange_1a) = zcolrate_ocnv(start_subrange_1a) + zj3n3(2) / pc_ocnv
4494!
4495!--       Total sink for organic vapor
4496          zcs_ocnv = zcs_tot + zj3n3(2) / pc_ocnv
4497!
4498!--       New gas phase concentration (#/m3)
4499          zcvap_new2 = pc_ocnv / ( 1.0_wp + ptstep * zcs_ocnv )
4500!
4501!--       Change in gas concentration (#/m3)
4502          zdvap2 = pc_ocnv - zcvap_new2
4503!
4504!--       Updated vapour concentration (#/m3)
4505          pc_ocnv = zcvap_new2
4506!
4507!--       Volume change of particles (m3(OC)/m3(air))
4508          zdvoloc = zcolrate_ocnv(start_subrange_1a:end_subrange_2b) / zcs_ocnv * amvoc * zdvap2
4509!
4510!--       Change of volume due to condensation in 1a-2b
4511          paero(start_subrange_1a:end_subrange_2b)%volc(2) =                                       &
4512                                          paero(start_subrange_1a:end_subrange_2b)%volc(2) + zdvoloc
4513!
4514!--       Change of number concentration in the smallest bin caused by nucleation (Jacobson (2005),
4515!--       eq. (16.75)). If zxocnv = 0, then the chosen nucleation mechanism doesn't take into
4516!--       account the non-volatile organic vapors and thus the paero doesn't have to be updated.
4517          IF ( zxocnv > 0.0_wp )  THEN
4518             paero(start_subrange_1a)%numc = paero(start_subrange_1a)%numc + zn_vs_c *             &
4519                                             zdvoloc(start_subrange_1a) / amvoc / ( n3 * zxocnv )
4520          ENDIF
4521       ENDIF
4522!
4523!--    5.1.2) Semivolatile organic compound: all bins except subrange 1
4524       zcs_ocsv = SUM( zcolrate(start_subrange_2a:end_subrange_2b) ) !< sink for semi-volatile organics
4525       IF ( pcocsv > 1.0E+10_wp  .AND.  zcs_ocsv > 1.0E-30  .AND. is_used( prtcl,'OC') )  THEN
4526!
4527!--       New gas phase concentration (#/m3)
4528          zcvap_new3 = pcocsv / ( 1.0_wp + ptstep * zcs_ocsv )
4529!
4530!--       Change in gas concentration (#/m3)
4531          zdvap3 = pcocsv - zcvap_new3
4532!
4533!--       Updated gas concentration (#/m3)
4534          pcocsv = zcvap_new3
4535!
4536!--       Volume change of particles (m3(OC)/m3(air))
4537          ss = start_subrange_2a
4538          ee = end_subrange_2b
4539          zdvoloc(ss:ee) = zdvoloc(ss:ee) + zcolrate(ss:ee) / zcs_ocsv * amvoc * zdvap3
4540!
4541!--       Change of volume due to condensation in 1a-2b
4542          paero(start_subrange_1a:end_subrange_2b)%volc(2) =                                       &
4543                                          paero(start_subrange_1a:end_subrange_2b)%volc(2) + zdvoloc
4544       ENDIF
4545!
4546!--    5.2) Sulphate: condensed on all bins
4547       IF ( pc_sa > 1.0E+10_wp  .AND.  zcs_tot > 1.0E-30_wp  .AND.  index_so4 > 0 )  THEN
4548!
4549!--    Ratio of mass transfer between nucleation and condensation
4550          zn_vs_c = 0.0_wp
4551          IF ( zj3n3(1) > 1.0_wp )  THEN
4552             zn_vs_c = ( zj3n3(1) ) / ( zj3n3(1) + pc_sa * zcolrate(start_subrange_1a) )
4553          ENDIF
4554!
4555!--       Collision rate in the smallest bin, including nucleation and condensation (see
4556!--       Jacobson (2005), eq. (16.73))
4557          zcolrate(start_subrange_1a) = zcolrate(start_subrange_1a) + zj3n3(1) / pc_sa
4558!
4559!--       Total sink for sulfate (1/s)
4560          zcs_su = zcs_tot + zj3n3(1) / pc_sa
4561!
4562!--       Sulphuric acid:
4563!--       New gas phase concentration (#/m3)
4564          zcvap_new1 = pc_sa / ( 1.0_wp + ptstep * zcs_su )
4565!
4566!--       Change in gas concentration (#/m3)
4567          zdvap1 = pc_sa - zcvap_new1
4568!
4569!--       Updating vapour concentration (#/m3)
4570          pc_sa = zcvap_new1
4571!
4572!--       Volume change of particles (m3(SO4)/m3(air)) by condensation
4573          zdvolsa = zcolrate(start_subrange_1a:end_subrange_2b) / zcs_su * amvh2so4 * zdvap1
4574!
4575!--       Change of volume concentration of sulphate in aerosol [fxm]
4576          paero(start_subrange_1a:end_subrange_2b)%volc(1) =                                       &
4577                                          paero(start_subrange_1a:end_subrange_2b)%volc(1) + zdvolsa
4578!
4579!--       Change of number concentration in the smallest bin caused by nucleation
4580!--       (Jacobson (2005), equation (16.75))
4581          IF ( zxsa > 0.0_wp )  THEN
4582             paero(start_subrange_1a)%numc = paero(start_subrange_1a)%numc + zn_vs_c *             &
4583                                             zdvolsa(start_subrange_1a) / amvh2so4 / ( n3 * zxsa)
4584          ENDIF
4585       ENDIF
4586!
4587!--    Partitioning of H2O, HNO3, and NH3: Dissolutional growth
4588       IF ( lspartition  .AND.  ( pchno3 > 1.0E+10_wp  .OR.  pc_nh3 > 1.0E+10_wp ) )  THEN
4589          CALL gpparthno3( ppres, ptemp, paero, pchno3, pc_nh3, pcw, pcs, zbeta, ptstep )
4590       ENDIF
4591    ENDIF
4592!
4593!-- Condensation of water vapour
4594    IF ( lscndh2oae )  THEN
4595       CALL gpparth2o( paero, ptemp, ppres, pcs, pcw, ptstep )
4596    ENDIF
4597
4598 END SUBROUTINE condensation
4599
4600!------------------------------------------------------------------------------!
4601! Description:
4602! ------------
4603!> Calculates the particle number and volume increase, and gas-phase
4604!> concentration decrease due to nucleation subsequent growth to detectable size
4605!> of 3 nm.
4606!
4607!> Method:
4608!> When the formed clusters grow by condensation (possibly also by self-
4609!> coagulation), their number is reduced due to scavenging to pre-existing
4610!> particles. Thus, the apparent nucleation rate at 3 nm is significantly lower
4611!> than the real nucleation rate (at ~1 nm).
4612!
4613!> Calculation of the formation rate of detectable particles at 3 nm (i.e. J3):
4614!> nj3 = 1: Kerminen, V.-M. and Kulmala, M. (2002), J. Aerosol Sci.,33, 609-622.
4615!> nj3 = 2: Lehtinen et al. (2007), J. Aerosol Sci., 38(9), 988-994.
4616!> nj3 = 3: Anttila et al. (2010), J. Aerosol Sci., 41(7), 621-636.
4617!
4618!> c = aerosol of critical radius (1 nm)
4619!> x = aerosol with radius 3 nm
4620!> 2 = wet or mean droplet
4621!
4622!> Called from subroutine condensation (in module salsa_dynamics_mod.f90)
4623!
4624!> Calls one of the following subroutines:
4625!>  - binnucl
4626!>  - ternucl
4627!>  - kinnucl
4628!>  - actnucl
4629!
4630!> fxm: currently only sulphuric acid grows particles from 1 to 3 nm
4631!>  (if asked from Markku, this is terribly wrong!!!)
4632!
4633!> Coded by:
4634!> Hannele Korhonen (FMI) 2005
4635!> Harri Kokkola (FMI) 2006
4636!> Matti Niskanen(FMI) 2012
4637!> Anton Laakso  (FMI) 2013
4638!------------------------------------------------------------------------------!
4639
4640 SUBROUTINE nucleation( paero, ptemp, prh, ppres, pc_sa, pc_ocnv, pc_nh3, ptstep, pj3n3, pxsa,     &
4641                        pxocnv )
4642
4643    IMPLICIT NONE
4644
4645    INTEGER(iwp) ::  iteration
4646
4647    REAL(wp) ::  zc_h2so4     !< H2SO4 conc. (#/cm3) !UNITS!
4648    REAL(wp) ::  zc_org       !< organic vapour conc. (#/cm3)
4649    REAL(wp) ::  zcc_c        !< Cunningham correct factor for c = critical (1nm)
4650    REAL(wp) ::  zcc_x        !< Cunningham correct factor for x = 3nm
4651    REAL(wp) ::  zcoags_c     !< coagulation sink (1/s) for c = critical (1nm)
4652    REAL(wp) ::  zcoags_x     !< coagulation sink (1/s) for x = 3nm
4653    REAL(wp) ::  zcoagstot    !< total particle losses due to coagulation, including condensation
4654                              !< and self-coagulation
4655    REAL(wp) ::  zcocnv_local !< organic vapour conc. (#/m3)
4656    REAL(wp) ::  zcsink       !< condensational sink (#/m2)
4657    REAL(wp) ::  zcsa_local   !< H2SO4 conc. (#/m3)
4658    REAL(wp) ::  zcv_c        !< mean relative thermal velocity (m/s) for c = critical (1nm)
4659    REAL(wp) ::  zcv_x        !< mean relative thermal velocity (m/s) for x = 3nm
4660    REAL(wp) ::  zdcrit       !< diameter of critical cluster (m)
4661    REAL(wp) ::  zdelta_vap   !< change of H2SO4 and organic vapour concentration (#/m3)
4662    REAL(wp) ::  zdfvap       !< air diffusion coefficient (m2/s)
4663    REAL(wp) ::  zdmean       !< mean diameter of existing particles (m)
4664    REAL(wp) ::  zeta         !< constant: proportional to ratio of CS/GR (m)
4665                              !< (condensation sink / growth rate)
4666    REAL(wp) ::  zgamma       !< proportionality factor ((nm2*m2)/h)
4667    REAL(wp) ::  z_gr_clust   !< growth rate of formed clusters (nm/h)
4668    REAL(wp) ::  z_gr_tot     !< total growth rate
4669    REAL(wp) ::  zj3          !< number conc. of formed 3nm particles (#/m3)
4670    REAL(wp) ::  zjnuc        !< nucleation rate at ~1nm (#/m3s)
4671    REAL(wp) ::  z_k_eff      !< effective cogulation coefficient for freshly nucleated particles
4672    REAL(wp) ::  zknud_c      !< Knudsen number for c = critical (1nm)
4673    REAL(wp) ::  zknud_x      !< Knudsen number for x = 3nm
4674    REAL(wp) ::  zkocnv       !< lever: zkocnv=1 --> organic compounds involved in nucleation
4675    REAL(wp) ::  zksa         !< lever: zksa=1 --> H2SO4 involved in nucleation
4676    REAL(wp) ::  zlambda      !< parameter for adjusting the growth rate due to self-coagulation
4677    REAL(wp) ::  zm_c         !< particle mass (kg) for c = critical (1nm)
4678    REAL(wp) ::  zm_para      !< Parameter m for calculating the coagulation sink (Eq. 5&6 in
4679                              !< Lehtinen et al. 2007)
4680    REAL(wp) ::  zm_x         !< particle mass (kg) for x = 3nm
4681    REAL(wp) ::  zmfp         !< mean free path of condesing vapour(m)
4682    REAL(wp) ::  zmixnh3      !< ammonia mixing ratio (ppt)
4683    REAL(wp) ::  zmyy         !< gas dynamic viscosity (N*s/m2)
4684    REAL(wp) ::  z_n_nuc      !< number of clusters/particles at the size range d1-dx (#/m3)
4685    REAL(wp) ::  znoc         !< number of organic molecules in critical cluster
4686    REAL(wp) ::  znsa         !< number of H2SO4 molecules in critical cluster
4687
4688    REAL(wp), INTENT(in) ::  pc_nh3   !< ammonia concentration (#/m3)
4689    REAL(wp), INTENT(in) ::  pc_ocnv  !< conc. of non-volatile OC (#/m3)
4690    REAL(wp), INTENT(in) ::  pc_sa    !< sulphuric acid conc. (#/m3)
4691    REAL(wp), INTENT(in) ::  ppres    !< ambient air pressure (Pa)
4692    REAL(wp), INTENT(in) ::  prh      !< ambient rel. humidity [0-1]
4693    REAL(wp), INTENT(in) ::  ptemp    !< ambient temperature (K)
4694    REAL(wp), INTENT(in) ::  ptstep   !< time step (s) of SALSA
4695
4696    REAL(wp), INTENT(inout) ::  pj3n3(2) !< formation mass rate of molecules (molec/m3s) for
4697                                         !< 1: H2SO4 and 2: organic vapour
4698
4699    REAL(wp), INTENT(out) ::  pxocnv  !< ratio of non-volatile organic vapours in 3 nm particles
4700    REAL(wp), INTENT(out) ::  pxsa    !< ratio of H2SO4 in 3 nm aerosol particles
4701
4702    REAL(wp), DIMENSION(nbins_aerosol) ::  zbeta       !< transitional correction factor
4703    REAL(wp), DIMENSION(nbins_aerosol) ::  zcc_2       !< Cunningham correct factor:2
4704    REAL(wp), DIMENSION(nbins_aerosol) ::  zcv_2       !< mean relative thermal velocity (m/s): 2
4705    REAL(wp), DIMENSION(nbins_aerosol) ::  zcv_c2      !< average velocity after coagulation: c & 2
4706    REAL(wp), DIMENSION(nbins_aerosol) ::  zcv_x2      !< average velocity after coagulation: x & 2
4707    REAL(wp), DIMENSION(nbins_aerosol) ::  zdc_2       !< particle diffusion coefficient (m2/s): 2
4708    REAL(wp), DIMENSION(nbins_aerosol) ::  zdc_c       !< particle diffusion coefficient (m2/s): c
4709    REAL(wp), DIMENSION(nbins_aerosol) ::  zdc_c2      !< sum of diffusion coef. for c and 2
4710    REAL(wp), DIMENSION(nbins_aerosol) ::  zdc_x       !< particle diffusion coefficient (m2/s): x
4711    REAL(wp), DIMENSION(nbins_aerosol) ::  zdc_x2      !< sum of diffusion coef. for: x & 2
4712    REAL(wp), DIMENSION(nbins_aerosol) ::  zgamma_f_2  !< zgamma_f for calculating zomega
4713    REAL(wp), DIMENSION(nbins_aerosol) ::  zgamma_f_c  !< zgamma_f for calculating zomega
4714    REAL(wp), DIMENSION(nbins_aerosol) ::  zgamma_f_x  !< zgamma_f for calculating zomega
4715    REAL(wp), DIMENSION(nbins_aerosol) ::  z_k_c2      !< coagulation coef. in the continuum
4716                                                       !< regime: c & 2
4717    REAL(wp), DIMENSION(nbins_aerosol) ::  z_k_x2      !< coagulation coef. in the continuum
4718                                                       !< regime: x & 2
4719    REAL(wp), DIMENSION(nbins_aerosol) ::  zknud       !< particle Knudsen number
4720    REAL(wp), DIMENSION(nbins_aerosol) ::  zknud_2     !< particle Knudsen number: 2
4721    REAL(wp), DIMENSION(nbins_aerosol) ::  zm_2        !< particle mass (kg): 2
4722    REAL(wp), DIMENSION(nbins_aerosol) ::  zomega_2c   !< zomega (m) for calculating zsigma: c & 2
4723    REAL(wp), DIMENSION(nbins_aerosol) ::  zomega_2x   !< zomega (m) for calculating zsigma: x & 2
4724    REAL(wp), DIMENSION(nbins_aerosol) ::  zomega_c    !< zomega (m) for calculating zsigma: c
4725    REAL(wp), DIMENSION(nbins_aerosol) ::  zomega_x    !< zomega (m) for calculating zsigma: x
4726    REAL(wp), DIMENSION(nbins_aerosol) ::  z_r_c2      !< sum of the radii: c & 2
4727    REAL(wp), DIMENSION(nbins_aerosol) ::  z_r_x2      !< sum of the radii: x & 2
4728    REAL(wp), DIMENSION(nbins_aerosol) ::  zsigma_c2   !<
4729    REAL(wp), DIMENSION(nbins_aerosol) ::  zsigma_x2   !<
4730
4731    TYPE(t