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

Last change on this file since 4598 was 4535, checked in by raasch, 11 months ago

bugfix for restart data format query

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