Home

Context Navigation

salsa_mod.f90 @ 3636

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