Home

Context Navigation

salsa_mod.f90 @ 4102

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