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

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

salsa_mod: time index error in salsa_emission_setup, lod=2

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