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

Last change on this file since 4012 was 4012, checked in by monakurppa, 2 years ago

remove salsa_util_mod.f90 and correct some bugs in salsa and salsa test case

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