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

Last change on this file since 4069 was 4069, checked in by Giersch, 2 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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