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

Last change on this file since 4286 was 4286, checked in by resler, 20 months ago

Merge branch resler into trunk

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