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

Last change on this file since 4131 was 4131, checked in by monakurppa, 23 months ago

Several changes in the salsa aerosol module:

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