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

Last change on this file since 4227 was 4227, checked in by gronemeier, 22 months ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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