Home

Context Navigation

salsa_mod.f90 @ 3483

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