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

Last change on this file since 4270 was 4270, checked in by monakurppa, 20 months ago

Updates for the aerosol module salsa

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