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

Last change on this file since 4525 was 4525, checked in by raasch, 8 months ago

salsa: added reading/writing of global restart data + reading/writing restart data with MPI-IO, bugfix for MPI-IO in plant_canopy_model_mod, tutorial user-interface dispersion_eularian_and_lpm_extended updated

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