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

Last change on this file since 4102 was 4102, checked in by suehring, 2 years ago

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

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