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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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