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

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