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

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

Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species

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