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

Last change on this file since 4527 was 4527, checked in by monakurppa, 8 months ago

salsa_mod: bug fixes in salsa_wrd_global, salsa_check_data_output and TEST/cases/urban_environment_salsa

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