Home

Context Navigation

salsa_mod.f90 @ 4118

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