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

Last change on this file since 4255 was 4255, checked in by monakurppa, 21 months ago

Bug corrected in salsa_mod in writing the run control

  • Property svn:keywords set to Id
File size: 534.8 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 4255 2019-10-04 11:50:55Z monakurppa $
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 )  nz, ny, nx, ncomponents_mass*nbins_aerosol, advect_particle_water
1096    ELSE
1097       WRITE( io, 5 )  nz, ny, nx, ncc*nbins_aerosol, advect_particle_water
1098    ENDIF
1099    IF ( .NOT. salsa_gases_from_chem )  THEN
1100       WRITE( io, 6 )  nz, ny, nx, ngases_salsa, salsa_gases_from_chem
1101    ENDIF
1102    WRITE( io, 7 )
1103    IF ( nsnucl > 0 )   WRITE( io, 8 ) nsnucl, nj3
1104    IF ( nlcoag )       WRITE( io, 9 )
1105    IF ( nlcnd )        WRITE( io, 10 ) nlcndgas, nlcndh2oae
1106    IF ( lspartition )  WRITE( io, 11 )
1107    IF ( nldepo )       WRITE( io, 12 ) nldepo_pcm, nldepo_surf
1108    WRITE( io, 13 )  reglim, nbin, bin_low_limits
1109    IF ( init_aerosol_type == 0 )  WRITE( io, 14 ) nsect
1110    WRITE( io, 15 ) ncc, listspec, mass_fracs_a, mass_fracs_b
1111    IF ( .NOT. salsa_gases_from_chem )  THEN
1112       WRITE( io, 16 ) ngases_salsa, h2so4_init, hno3_init, nh3_init, ocnv_init, ocsv_init
1113    ENDIF
1114    WRITE( io, 17 )  init_aerosol_type, init_gases_type
1115    IF ( init_aerosol_type == 0 )  THEN
1116       WRITE( io, 18 )  dpg, sigmag, n_lognorm
1117    ELSE
1118       WRITE( io, 19 )
1119    ENDIF
1120    IF ( nest_salsa )  WRITE( io, 20 )  nest_salsa
1121    WRITE( io, 21 ) salsa_emission_mode
1122    IF ( salsa_emission_mode == 'uniform' )  THEN
1123       WRITE( io, 22 ) surface_aerosol_flux, aerosol_flux_dpg, aerosol_flux_sigmag,                &
1124                       aerosol_flux_mass_fracs_a
1125    ENDIF
1126    IF ( SUM( aerosol_flux_mass_fracs_b ) > 0.0_wp  .OR. salsa_emission_mode == 'read_from_file' ) &
1127    THEN
1128       WRITE( io, 23 )
1129    ENDIF
1130
11311   FORMAT (//' SALSA information:'/                                                               &
1132              ' ------------------------------'/)
11332   FORMAT   ('    Starts at: skip_time_do_salsa = ', F10.2, '  s')
11343   FORMAT  (/'    Timestep: dt_salsa = ', F6.2, '  s')
11354   FORMAT  (/'    Array shape (z,y,x,bins):'/                                                     &
1136              '       aerosol_number:  ', 4(I3)) 
11375   FORMAT  (/'       aerosol_mass:    ', 4(I3),/                                                  &
1138              '       (advect_particle_water = ', L1, ')')
11396   FORMAT   ('       salsa_gas: ', 4(I3),/                                                        &
1140              '       (salsa_gases_from_chem = ', L1, ')')
11417   FORMAT  (/'    Aerosol dynamic processes included: ')
11428   FORMAT  (/'       nucleation (scheme = ', I1, ' and J3 parametrization = ', I1, ')')
11439   FORMAT  (/'       coagulation')
114410  FORMAT  (/'       condensation (of precursor gases = ', L1, ' and water vapour = ', L1, ')' )
114511  FORMAT  (/'       dissolutional growth by HNO3 and NH3')
114612  FORMAT  (/'       dry deposition (on vegetation = ', L1, ' and on topography = ', L1, ')')
114713  FORMAT  (/'    Aerosol bin subrange limits (in metres): ',  3(ES10.2E3), /                     &
1148              '    Number of size bins for each aerosol subrange: ', 2I3,/                         &
1149              '    Aerosol bin limits (in metres): ', 9(ES10.2E3))
115014  FORMAT   ('    Initial number concentration in bins at the lowest level (#/m**3):', 9(ES10.2E3))
115115  FORMAT  (/'    Number of chemical components used: ', I1,/                                     &
1152              '       Species: ',7(A6),/                                                           &
1153              '    Initial relative contribution of each species to particle volume in:',/         &
1154              '       a-bins: ', 7(F6.3),/                                                         &
1155              '       b-bins: ', 7(F6.3))
115616  FORMAT  (/'    Number of gaseous tracers used: ', I1,/                                         &
1157              '    Initial gas concentrations:',/                                                  &
1158              '       H2SO4: ',ES12.4E3, ' #/m**3',/                                               &
1159              '       HNO3:  ',ES12.4E3, ' #/m**3',/                                               &
1160              '       NH3:   ',ES12.4E3, ' #/m**3',/                                               &
1161              '       OCNV:  ',ES12.4E3, ' #/m**3',/                                               &
1162              '       OCSV:  ',ES12.4E3, ' #/m**3')
116317   FORMAT (/'   Initialising concentrations: ', /                                                &
1164              '      Aerosol size distribution: init_aerosol_type = ', I1,/                        &
1165              '      Gas concentrations: init_gases_type = ', I1 )
116618   FORMAT ( '      Mode diametres: dpg(nmod) = ', 7(F7.3), ' (m)', /                             &
1167              '      Standard deviation: sigmag(nmod) = ', 7(F7.2),/                               &
1168              '      Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3), ' (#/m3)' )
116919   FORMAT (/'      Size distribution read from a file.')
117020   FORMAT (/'   Nesting for salsa variables: ', L1 )
117121   FORMAT (/'   Emissions: salsa_emission_mode = ', A )
117222   FORMAT (/'      surface_aerosol_flux = ', ES12.4E3, ' #/m**2/s', /                            &
1173              '      aerosol_flux_dpg     =  ', 7(F7.3), ' (m)', /                                 &
1174              '      aerosol_flux_sigmag  =  ', 7(F7.2), /                                         &
1175              '      aerosol_mass_fracs_a =  ', 7(ES12.4E3) )
117623   FORMAT (/'      (currently all emissions are soluble!)')
1177
1178 END SUBROUTINE salsa_header
1179
1180!------------------------------------------------------------------------------!
1181! Description:
1182! ------------
1183!> Allocate SALSA arrays and define pointers if required
1184!------------------------------------------------------------------------------!
1185 SUBROUTINE salsa_init_arrays
1186
1187    USE advec_ws,                                                                                  &
1188        ONLY: ws_init_flags_scalar
1189
1190    USE surface_mod,                                                                               &
1191        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
1192
1193    IMPLICIT NONE
1194
1195    INTEGER(iwp) ::  gases_available !< Number of available gas components in the chemistry model
1196    INTEGER(iwp) ::  i               !< loop index for allocating
1197    INTEGER(iwp) ::  ii              !< index for indexing chemical components
1198    INTEGER(iwp) ::  l               !< loop index for allocating: surfaces
1199    INTEGER(iwp) ::  lsp             !< loop index for chem species in the chemistry model
1200
1201    gases_available = 0
1202!
1203!-- Allocate prognostic variables (see salsa_swap_timelevel)
1204!
1205!-- Set derived indices:
1206!-- (This does the same as the subroutine salsa_initialize in SALSA/UCLALES-SALSA)
1207    start_subrange_1a = 1  ! 1st index of subrange 1a
1208    start_subrange_2a = start_subrange_1a + nbin(1)  ! 1st index of subrange 2a
1209    end_subrange_1a   = start_subrange_2a - 1        ! last index of subrange 1a
1210    end_subrange_2a   = end_subrange_1a + nbin(2)    ! last index of subrange 2a
1211
1212!
1213!-- If the fraction of insoluble aerosols in subrange 2 is zero: do not allocate arrays for them
1214    IF ( nf2a > 0.999999_wp  .AND.  SUM( mass_fracs_b ) < 0.00001_wp )  THEN
1215       no_insoluble = .TRUE.
1216       start_subrange_2b = end_subrange_2a+1  ! 1st index of subrange 2b
1217       end_subrange_2b   = end_subrange_2a    ! last index of subrange 2b
1218    ELSE
1219       start_subrange_2b = start_subrange_2a + nbin(2)  ! 1st index of subrange 2b
1220       end_subrange_2b   = end_subrange_2a + nbin(2)    ! last index of subrange 2b
1221    ENDIF
1222
1223    nbins_aerosol = end_subrange_2b   ! total number of aerosol size bins
1224!
1225!-- Create index tables for different aerosol components
1226    CALL component_index_constructor( prtcl, ncc, maxspec, listspec )
1227
1228    ncomponents_mass = ncc
1229    IF ( advect_particle_water )  ncomponents_mass = ncc + 1  ! Add water
1230!
1231!-- Indices for chemical components used (-1 = not used)
1232    ii = 0
1233    IF ( is_used( prtcl, 'SO4' ) )  THEN
1234       index_so4 = get_index( prtcl,'SO4' )
1235       ii = ii + 1
1236    ENDIF
1237    IF ( is_used( prtcl,'OC' ) )  THEN
1238       index_oc = get_index(prtcl, 'OC')
1239       ii = ii + 1
1240    ENDIF
1241    IF ( is_used( prtcl, 'BC' ) )  THEN
1242       index_bc = get_index( prtcl, 'BC' )
1243       ii = ii + 1
1244    ENDIF
1245    IF ( is_used( prtcl, 'DU' ) )  THEN
1246       index_du = get_index( prtcl, 'DU' )
1247       ii = ii + 1
1248    ENDIF
1249    IF ( is_used( prtcl, 'SS' ) )  THEN
1250       index_ss = get_index( prtcl, 'SS' )
1251       ii = ii + 1
1252    ENDIF
1253    IF ( is_used( prtcl, 'NO' ) )  THEN
1254       index_no = get_index( prtcl, 'NO' )
1255       ii = ii + 1
1256    ENDIF
1257    IF ( is_used( prtcl, 'NH' ) )  THEN
1258       index_nh = get_index( prtcl, 'NH' )
1259       ii = ii + 1
1260    ENDIF
1261!
1262!-- All species must be known
1263    IF ( ii /= ncc )  THEN
1264       message_string = 'Unknown aerosol species/component(s) given in the initialization'
1265       CALL message( 'salsa_mod: salsa_init', 'PA0600', 1, 2, 0, 6, 0 )
1266    ENDIF
1267!
1268!-- Allocate:
1269    ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass),                    &
1270              bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),&
1271              nsect(nbins_aerosol), massacc(nbins_aerosol) )
1272    ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) )
1273    IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1274    ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1275!
1276!-- Initialise the sectional particle size distribution
1277    CALL set_sizebins
1278!
1279!-- Aerosol number concentration
1280    ALLOCATE( aerosol_number(nbins_aerosol) )
1281    ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
1282              nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
1283              nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1284    nconc_1 = 0.0_wp
1285    nconc_2 = 0.0_wp
1286    nconc_3 = 0.0_wp
1287
1288    DO i = 1, nbins_aerosol
1289       aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => nconc_1(:,:,:,i)
1290       aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => nconc_2(:,:,:,i)
1291       aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i)
1292       ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                         &
1293                 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                         &
1294                 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
1295                 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
1296                 aerosol_number(i)%init(nzb:nzt+1),                                                &
1297                 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1298       aerosol_number(i)%init = nclim
1299       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
1300          ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
1301          aerosol_number(i)%source = 0.0_wp
1302       ENDIF
1303    ENDDO
1304
1305!
1306!-- Aerosol mass concentration
1307    ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) )
1308    ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
1309              mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
1310              mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) )
1311    mconc_1 = 0.0_wp
1312    mconc_2 = 0.0_wp
1313    mconc_3 = 0.0_wp
1314
1315    DO i = 1, ncomponents_mass*nbins_aerosol
1316       aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => mconc_1(:,:,:,i)
1317       aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => mconc_2(:,:,:,i)
1318       aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i)
1319       ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                           &
1320                 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                           &
1321                 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
1322                 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
1323                 aerosol_mass(i)%init(nzb:nzt+1),                                                  &
1324                 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
1325       aerosol_mass(i)%init = mclim
1326       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
1327          ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
1328          aerosol_mass(i)%source = 0.0_wp
1329       ENDIF
1330    ENDDO
1331
1332!
1333!-- Surface fluxes: answs = aerosol number, amsws = aerosol mass
1334!
1335!-- Horizontal surfaces: default type
1336    DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
1337       ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) )
1338       ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) )
1339       surf_def_h(l)%answs = 0.0_wp
1340       surf_def_h(l)%amsws = 0.0_wp
1341    ENDDO
1342!
1343!-- Horizontal surfaces: natural type
1344    ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) )
1345    ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) )
1346    surf_lsm_h%answs = 0.0_wp
1347    surf_lsm_h%amsws = 0.0_wp
1348!
1349!-- Horizontal surfaces: urban type
1350    ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) )
1351    ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) )
1352    surf_usm_h%answs = 0.0_wp
1353    surf_usm_h%amsws = 0.0_wp
1354
1355!
1356!-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing
1357    DO  l = 0, 3
1358       ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) )
1359       surf_def_v(l)%answs = 0.0_wp
1360       ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1361       surf_def_v(l)%amsws = 0.0_wp
1362
1363       ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) )
1364       surf_lsm_v(l)%answs = 0.0_wp
1365       ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1366       surf_lsm_v(l)%amsws = 0.0_wp
1367
1368       ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) )
1369       surf_usm_v(l)%answs = 0.0_wp
1370       ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1371       surf_usm_v(l)%amsws = 0.0_wp
1372
1373    ENDDO
1374
1375!
1376!-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV)
1377!-- (number concentration (#/m3) )
1378!
1379!-- If chemistry is on, read gas phase concentrations from there. Otherwise,
1380!-- allocate salsa_gas array.
1381
1382    IF ( air_chemistry )  THEN
1383       DO  lsp = 1, nvar
1384          SELECT CASE ( TRIM( chem_species(lsp)%name ) )
1385             CASE ( 'H2SO4', 'h2so4' )
1386                gases_available = gases_available + 1
1387                gas_index_chem(1) = lsp
1388             CASE ( 'HNO3', 'hno3' )
1389                gases_available = gases_available + 1
1390                gas_index_chem(2) = lsp
1391             CASE ( 'NH3', 'nh3' )
1392                gases_available = gases_available + 1
1393                gas_index_chem(3) = lsp
1394             CASE ( 'OCNV', 'ocnv' )
1395                gases_available = gases_available + 1
1396                gas_index_chem(4) = lsp
1397             CASE ( 'OCSV', 'ocsv' )
1398                gases_available = gases_available + 1
1399                gas_index_chem(5) = lsp
1400          END SELECT
1401       ENDDO
1402
1403       IF ( gases_available == ngases_salsa )  THEN
1404          salsa_gases_from_chem = .TRUE.
1405       ELSE
1406          WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// &
1407                                     'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)'
1408       CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 )
1409       ENDIF
1410
1411    ELSE
1412
1413       ALLOCATE( salsa_gas(ngases_salsa) )
1414       ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
1415                 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
1416                 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
1417       gconc_1 = 0.0_wp
1418       gconc_2 = 0.0_wp
1419       gconc_3 = 0.0_wp
1420
1421       DO i = 1, ngases_salsa
1422          salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => gconc_1(:,:,:,i)
1423          salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => gconc_2(:,:,:,i)
1424          salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i)
1425          ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
1426                    salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
1427                    salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1428                    salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1429                    salsa_gas(i)%init(nzb:nzt+1),                              &
1430                    salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1431          salsa_gas(i)%init = nclim
1432          IF ( include_emission )  THEN
1433             ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
1434             salsa_gas(i)%source = 0.0_wp
1435          ENDIF
1436       ENDDO
1437!
1438!--    Surface fluxes: gtsws = gaseous tracer flux
1439!
1440!--    Horizontal surfaces: default type
1441       DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
1442          ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) )
1443          surf_def_h(l)%gtsws = 0.0_wp
1444       ENDDO
1445!--    Horizontal surfaces: natural type
1446       ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) )
1447       surf_lsm_h%gtsws = 0.0_wp
1448!--    Horizontal surfaces: urban type
1449       ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) )
1450       surf_usm_h%gtsws = 0.0_wp
1451!
1452!--    Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
1453!--    westward (l=3) facing
1454       DO  l = 0, 3
1455          ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) )
1456          surf_def_v(l)%gtsws = 0.0_wp
1457          ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) )
1458          surf_lsm_v(l)%gtsws = 0.0_wp
1459          ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) )
1460          surf_usm_v(l)%gtsws = 0.0_wp
1461       ENDDO
1462    ENDIF
1463
1464    IF ( ws_scheme_sca )  THEN
1465
1466       IF ( salsa )  THEN
1467          ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1468          sums_salsa_ws_l = 0.0_wp
1469       ENDIF
1470
1471    ENDIF
1472!
1473!-- Set control flags for decycling only at lateral boundary cores. Within the inner cores the
1474!-- decycle flag is set to .FALSE.. Even though it does not affect the setting of chemistry boundary
1475!-- conditions, this flag is used to set advection control flags appropriately.
1476    decycle_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0  .OR.  nxr == nx )
1477    decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0  .OR.  nyn == ny )
1478!
1479!-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have
1480!-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce
1481!-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
1482!-- applied near these boundaries. To get rid-off this, set-up additional flags that control the
1483!-- order of the scalar advection scheme near the lateral boundaries for passive scalars with
1484!-- decycling.
1485    IF ( scalar_advec == 'ws-scheme' )  THEN
1486       ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1487!
1488!--    In case of decycling, set Neuman boundary conditions for wall_flags_0 bit 31 instead of
1489!--    cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see
1490!--    the following comment). Note, since several also other modules may access this bit but may
1491!--    have other boundary conditions, the original value of wall_flags_0 bit 31 must not be
1492!--    modified. Hence, store the boundary conditions directly on salsa_advc_flags_s.
1493!--    salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used
1494!--    to control the numerical order.
1495!--    Initialize with flag 31 only.
1496       salsa_advc_flags_s = 0
1497       salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_0, 31 ) )
1498
1499       IF ( decycle_salsa_ns )  THEN
1500          IF ( nys == 0 )  THEN
1501             DO  i = 1, nbgp
1502                salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ),   &
1503                                                       IBCLR( salsa_advc_flags_s(:,nys,:), 31 ),   &
1504                                                       BTEST( salsa_advc_flags_s(:,nys,:), 31 ) )
1505             ENDDO
1506          ENDIF
1507          IF ( nyn == ny )  THEN
1508             DO  i = 1, nbgp
1509                salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ),   &
1510                                                       IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ),   &
1511                                                       BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) )
1512             ENDDO
1513          ENDIF
1514       ENDIF
1515       IF ( decycle_salsa_lr )  THEN
1516          IF ( nxl == 0 )  THEN
1517             DO  i = 1, nbgp
1518                salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ),   &
1519                                                       IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ),   &
1520                                                       BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) )
1521             ENDDO
1522          ENDIF
1523          IF ( nxr == nx )  THEN
1524             DO  i = 1, nbgp
1525                salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ),   &
1526                                                       IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ),   &
1527                                                       BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) )
1528             ENDDO
1529          ENDIF
1530       ENDIF
1531!
1532!--    To initialise the advection flags appropriately, pass the boundary flags to
1533!--    ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive
1534!--    scalar is being treated and the horizontal advection terms are degraded already 2 grid points
1535!--    before the lateral boundary. Also, extended degradation zones are applied, where
1536!--    horizontal advection of scalars is discretised by the first-order scheme at all grid points
1537!--    in the vicinity of buildings (<= 3 grid points). Even though no building is within the
1538!--    numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order
1539!--    of the horizontal advection scheme is successively upgraded.
1540!--    These degradations of the advection scheme are done to avoid stationary numerical
1541!--    oscillations, which are responsible for high concentration maxima that may appear e.g. under
1542!--    shear-free stable conditions.
1543       CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_salsa_lr,    &
1544                                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_salsa_ns,    &
1545                                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_salsa_lr,    &
1546                                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_salsa_ns,    &
1547                                  salsa_advc_flags_s, .TRUE. )
1548    ENDIF
1549
1550
1551 END SUBROUTINE salsa_init_arrays
1552
1553!------------------------------------------------------------------------------!
1554! Description:
1555! ------------
1556!> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA.
1557!> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are
1558!> also merged here.
1559!------------------------------------------------------------------------------!
1560 SUBROUTINE salsa_init
1561
1562    IMPLICIT NONE
1563
1564    INTEGER(iwp) :: i   !<
1565    INTEGER(iwp) :: ib  !< loop index for aerosol number bins
1566    INTEGER(iwp) :: ic  !< loop index for aerosol mass bins
1567    INTEGER(iwp) :: ig  !< loop index for gases
1568    INTEGER(iwp) :: j   !<
1569
1570    IF ( debug_output )  CALL debug_message( 'salsa_init', 'start' )
1571
1572    bin_low_limits = 0.0_wp
1573    k_topo_top     = 0
1574    nsect          = 0.0_wp
1575    massacc        = 1.0_wp
1576!
1577!-- Initialise
1578    IF ( nldepo )  sedim_vd = 0.0_wp
1579
1580    IF ( .NOT. salsa_gases_from_chem )  THEN
1581       IF ( .NOT. read_restart_data_salsa )  THEN
1582          salsa_gas(1)%conc = h2so4_init
1583          salsa_gas(2)%conc = hno3_init
1584          salsa_gas(3)%conc = nh3_init
1585          salsa_gas(4)%conc = ocnv_init
1586          salsa_gas(5)%conc = ocsv_init
1587       ENDIF
1588       DO  ig = 1, ngases_salsa
1589          salsa_gas(ig)%conc_p    = 0.0_wp
1590          salsa_gas(ig)%tconc_m   = 0.0_wp
1591          salsa_gas(ig)%flux_s    = 0.0_wp
1592          salsa_gas(ig)%diss_s    = 0.0_wp
1593          salsa_gas(ig)%flux_l    = 0.0_wp
1594          salsa_gas(ig)%diss_l    = 0.0_wp
1595          salsa_gas(ig)%sums_ws_l = 0.0_wp
1596          salsa_gas(ig)%conc_p    = salsa_gas(ig)%conc
1597       ENDDO
1598!
1599!--    Set initial value for gas compound tracer
1600       salsa_gas(1)%init = h2so4_init
1601       salsa_gas(2)%init = hno3_init
1602       salsa_gas(3)%init = nh3_init
1603       salsa_gas(4)%init = ocnv_init
1604       salsa_gas(5)%init = ocsv_init
1605    ENDIF
1606!
1607!-- Aerosol radius in each bin: dry and wet (m)
1608    ra_dry = 1.0E-10_wp
1609!
1610!-- Initialise location-dependent aerosol size distributions and chemical compositions:
1611    CALL aerosol_init
1612
1613!-- Initalisation run of SALSA + calculate the vertical top index of the topography
1614    DO  i = nxl, nxr
1615       DO  j = nys, nyn
1616
1617          k_topo_top(j,i) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,j,i), 12 ) ), DIM = 1 ) - 1
1618
1619          CALL salsa_driver( i, j, 1 )
1620          CALL salsa_diagnostics( i, j )
1621       ENDDO
1622    ENDDO
1623
1624    DO  ib = 1, nbins_aerosol
1625       aerosol_number(ib)%conc_p    = aerosol_number(ib)%conc
1626       aerosol_number(ib)%tconc_m   = 0.0_wp
1627       aerosol_number(ib)%flux_s    = 0.0_wp
1628       aerosol_number(ib)%diss_s    = 0.0_wp
1629       aerosol_number(ib)%flux_l    = 0.0_wp
1630       aerosol_number(ib)%diss_l    = 0.0_wp
1631       aerosol_number(ib)%sums_ws_l = 0.0_wp
1632    ENDDO
1633    DO  ic = 1, ncomponents_mass*nbins_aerosol
1634       aerosol_mass(ic)%conc_p    = aerosol_mass(ic)%conc
1635       aerosol_mass(ic)%tconc_m   = 0.0_wp
1636       aerosol_mass(ic)%flux_s    = 0.0_wp
1637       aerosol_mass(ic)%diss_s    = 0.0_wp
1638       aerosol_mass(ic)%flux_l    = 0.0_wp
1639       aerosol_mass(ic)%diss_l    = 0.0_wp
1640       aerosol_mass(ic)%sums_ws_l = 0.0_wp
1641    ENDDO
1642!
1643!
1644!-- Initialise the deposition scheme and surface types
1645    IF ( nldepo )  CALL init_deposition
1646
1647    IF ( include_emission )  THEN
1648!
1649!--    Read in and initialize emissions
1650       CALL salsa_emission_setup( .TRUE. )
1651       IF ( .NOT. salsa_gases_from_chem  .AND.  salsa_emission_mode == 'read_from_file' )  THEN
1652          CALL salsa_gas_emission_setup( .TRUE. )
1653       ENDIF
1654    ENDIF
1655!
1656!-- Partition and dissolutional growth by gaseous HNO3 and NH3
1657    IF ( index_no > 0  .AND.  index_nh > 0  .AND.  index_so4 > 0 )  lspartition = .TRUE.
1658
1659    IF ( debug_output )  CALL debug_message( 'salsa_init', 'end' )
1660
1661 END SUBROUTINE salsa_init
1662
1663!------------------------------------------------------------------------------!
1664! Description:
1665! ------------
1666!> Initializes particle size distribution grid by calculating size bin limits
1667!> and mid-size for *dry* particles in each bin. Called from salsa_initialize
1668!> (only at the beginning of simulation).
1669!> Size distribution described using:
1670!>   1) moving center method (subranges 1 and 2)
1671!>      (Jacobson, Atmos. Env., 31, 131-144, 1997)
1672!>   2) fixed sectional method (subrange 3)
1673!> Size bins in each subrange are spaced logarithmically
1674!> based on given subrange size limits and bin number.
1675!
1676!> Mona changed 06/2017: Use geometric mean diameter to describe the mean
1677!> particle diameter in a size bin, not the arithmeric mean which clearly
1678!> overestimates the total particle volume concentration.
1679!
1680!> Coded by:
1681!> Hannele Korhonen (FMI) 2005
1682!> Harri Kokkola (FMI) 2006
1683!
1684!> Bug fixes for box model + updated for the new aerosol datatype:
1685!> Juha Tonttila (FMI) 2014
1686!------------------------------------------------------------------------------!
1687 SUBROUTINE set_sizebins
1688
1689    IMPLICIT NONE
1690
1691    INTEGER(iwp) ::  cc  !< running index
1692    INTEGER(iwp) ::  dd  !< running index
1693
1694    REAL(wp) ::  ratio_d  !< ratio of the upper and lower diameter of subranges
1695
1696    aero(:)%dwet     = 1.0E-10_wp
1697    aero(:)%veqh2o   = 1.0E-10_wp
1698    aero(:)%numc     = nclim
1699    aero(:)%core     = 1.0E-10_wp
1700    DO  cc = 1, maxspec+1    ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
1701       aero(:)%volc(cc) = 0.0_wp
1702    ENDDO
1703!
1704!-- vlolim&vhilim: min & max *dry* volumes [fxm]
1705!-- dmid: bin mid *dry* diameter (m)
1706!-- vratiolo&vratiohi: volume ratio between the center and low/high limit
1707!
1708!-- 1) Size subrange 1:
1709    ratio_d = reglim(2) / reglim(1)   ! section spacing (m)
1710    DO  cc = start_subrange_1a, end_subrange_1a
1711       aero(cc)%vlolim = api6 * ( reglim(1) * ratio_d**( REAL( cc-1 ) / nbin(1) ) )**3
1712       aero(cc)%vhilim = api6 * ( reglim(1) * ratio_d**( REAL( cc ) / nbin(1) ) )**3
1713       aero(cc)%dmid = SQRT( ( aero(cc)%vhilim / api6 )**0.33333333_wp *                           &
1714                             ( aero(cc)%vlolim / api6 )**0.33333333_wp )
1715       aero(cc)%vratiohi = aero(cc)%vhilim / ( api6 * aero(cc)%dmid**3 )
1716       aero(cc)%vratiolo = aero(cc)%vlolim / ( api6 * aero(cc)%dmid**3 )
1717    ENDDO
1718!
1719!-- 2) Size subrange 2:
1720!-- 2.1) Sub-subrange 2a: high hygroscopicity
1721    ratio_d = reglim(3) / reglim(2)   ! section spacing
1722    DO  dd = start_subrange_2a, end_subrange_2a
1723       cc = dd - start_subrange_2a
1724       aero(dd)%vlolim = api6 * ( reglim(2) * ratio_d**( REAL( cc ) / nbin(2) ) )**3
1725       aero(dd)%vhilim = api6 * ( reglim(2) * ratio_d**( REAL( cc+1 ) / nbin(2) ) )**3
1726       aero(dd)%dmid = SQRT( ( aero(dd)%vhilim / api6 )**0.33333333_wp *                           &
1727                             ( aero(dd)%vlolim / api6 )**0.33333333_wp )
1728       aero(dd)%vratiohi = aero(dd)%vhilim / ( api6 * aero(dd)%dmid**3 )
1729       aero(dd)%vratiolo = aero(dd)%vlolim / ( api6 * aero(dd)%dmid**3 )
1730    ENDDO
1731!
1732!-- 2.2) Sub-subrange 2b: low hygroscopicity
1733    IF ( .NOT. no_insoluble )  THEN
1734       aero(start_subrange_2b:end_subrange_2b)%vlolim   = aero(start_subrange_2a:end_subrange_2a)%vlolim
1735       aero(start_subrange_2b:end_subrange_2b)%vhilim   = aero(start_subrange_2a:end_subrange_2a)%vhilim
1736       aero(start_subrange_2b:end_subrange_2b)%dmid     = aero(start_subrange_2a:end_subrange_2a)%dmid
1737       aero(start_subrange_2b:end_subrange_2b)%vratiohi = aero(start_subrange_2a:end_subrange_2a)%vratiohi
1738       aero(start_subrange_2b:end_subrange_2b)%vratiolo = aero(start_subrange_2a:end_subrange_2a)%vratiolo
1739    ENDIF
1740!
1741!-- Initialize the wet diameter with the bin dry diameter to avoid numerical problems later
1742    aero(:)%dwet = aero(:)%dmid
1743!
1744!-- Save bin limits (lower diameter) to be delivered to PALM if needed
1745    DO cc = 1, nbins_aerosol
1746       bin_low_limits(cc) = ( aero(cc)%vlolim / api6 )**0.33333333_wp
1747    ENDDO
1748
1749 END SUBROUTINE set_sizebins
1750
1751!------------------------------------------------------------------------------!
1752! Description:
1753! ------------
1754!> Initilize altitude-dependent aerosol size distributions and compositions.
1755!>
1756!> Mona added 06/2017: Correct the number and mass concentrations by normalizing
1757!< by the given total number and mass concentration.
1758!>
1759!> Tomi Raatikainen, FMI, 29.2.2016
1760!------------------------------------------------------------------------------!
1761 SUBROUTINE aerosol_init
1762
1763    USE netcdf_data_input_mod,                                                                     &
1764        ONLY:  check_existence, close_input_file, get_dimension_length,                            &
1765               get_attribute, get_variable,                                                        &
1766               inquire_num_variables, inquire_variable_names,                                      &
1767               open_read_file
1768
1769    IMPLICIT NONE
1770
1771    CHARACTER(LEN=25),  DIMENSION(:), ALLOCATABLE ::  cc_name    !< chemical component name
1772    CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names
1773
1774    INTEGER(iwp) ::  ee        !< index: end
1775    INTEGER(iwp) ::  i         !< loop index: x-direction
1776    INTEGER(iwp) ::  ib        !< loop index: size bins
1777    INTEGER(iwp) ::  ic        !< loop index: chemical components
1778    INTEGER(iwp) ::  id_dyn    !< NetCDF id of PIDS_DYNAMIC_SALSA
1779    INTEGER(iwp) ::  ig        !< loop index: gases
1780    INTEGER(iwp) ::  j         !< loop index: y-direction
1781    INTEGER(iwp) ::  k         !< loop index: z-direction
1782    INTEGER(iwp) ::  lod_aero  !< level of detail of inital aerosol concentrations
1783    INTEGER(iwp) ::  num_vars  !< number of variables
1784    INTEGER(iwp) ::  pr_nbins  !< number of aerosol size bins in file
1785    INTEGER(iwp) ::  pr_ncc    !< number of aerosol chemical components in file
1786    INTEGER(iwp) ::  pr_nz     !< number of vertical grid-points in file
1787    INTEGER(iwp) ::  prunmode  !< running mode of SALSA
1788    INTEGER(iwp) ::  ss        !< index: start
1789
1790    INTEGER(iwp), DIMENSION(maxspec) ::  cc_in2mod
1791
1792    LOGICAL  ::  netcdf_extend = .FALSE. !< Flag: netcdf file exists
1793
1794    REAL(wp) ::  flag  !< flag to mask topography grid points
1795
1796    REAL(wp), DIMENSION(nbins_aerosol) ::  core   !< size of the bin mid aerosol particle
1797
1798    REAL(wp), DIMENSION(0:nz+1) ::  pnf2a   !< number fraction in 2a
1799    REAL(wp), DIMENSION(0:nz+1) ::  pmfoc1a !< mass fraction of OC in 1a
1800
1801    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol)   ::  pndist  !< vertical profile of size dist. (#/m3)
1802    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2a   !< mass distributions in subrange 2a
1803    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2b   !< mass distributions in subrange 2b
1804
1805    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_dmid  !< vertical profile of aerosol bin diameters
1806    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_z     !< z levels of profiles
1807
1808    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_a  !< mass fraction: a
1809    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_b  !< and b
1810
1811    cc_in2mod = 0
1812    prunmode = 1
1813!
1814!-- Bin mean aerosol particle volume (m3)
1815    core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3
1816!
1817!-- Set concentrations to zero
1818    pndist(:,:)  = 0.0_wp
1819    pnf2a(:)     = nf2a
1820    pmf2a(:,:)   = 0.0_wp
1821    pmf2b(:,:)   = 0.0_wp
1822    pmfoc1a(:)   = 0.0_wp
1823
1824    IF ( init_aerosol_type == 1 )  THEN
1825!
1826!--    Read input profiles from PIDS_DYNAMIC_SALSA
1827#if defined( __netcdf )
1828!
1829!--    Location-dependent size distributions and compositions.
1830       INQUIRE( FILE = TRIM( input_file_dynamic ) //  TRIM( coupling_char ), EXIST = netcdf_extend )
1831       IF ( netcdf_extend )  THEN
1832!
1833!--       Open file in read-only mode
1834          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
1835!
1836!--       At first, inquire all variable names
1837          CALL inquire_num_variables( id_dyn, num_vars )
1838!
1839!--       Allocate memory to store variable names
1840          ALLOCATE( var_names(1:num_vars) )
1841          CALL inquire_variable_names( id_dyn, var_names )
1842!
1843!--       Inquire vertical dimension and number of aerosol chemical components
1844          CALL get_dimension_length( id_dyn, pr_nz, 'z' )
1845          IF ( pr_nz /= nz )  THEN
1846             WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//&
1847                                        'the number of numeric grid points.'
1848             CALL message( 'aerosol_init', 'PA0601', 1, 2, 0, 6, 0 )
1849          ENDIF
1850          CALL get_dimension_length( id_dyn, pr_ncc, 'composition_index' )
1851!
1852!--       Allocate memory
1853          ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                              &
1854                    pr_mass_fracs_b(nzb:nzt+1,pr_ncc) )
1855          pr_mass_fracs_a = 0.0_wp
1856          pr_mass_fracs_b = 0.0_wp
1857!
1858!--       Read vertical levels
1859          CALL get_variable( id_dyn, 'z', pr_z )
1860!
1861!--       Read the names of chemical components
1862          IF ( check_existence( var_names, 'composition_name' ) )  THEN
1863             CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc )
1864          ELSE
1865             WRITE( message_string, * ) 'Missing composition_name in ' // TRIM( input_file_dynamic )
1866             CALL message( 'aerosol_init', 'PA0655', 1, 2, 0, 6, 0 )
1867          ENDIF
1868!
1869!--       Define the index of each chemical component in the model
1870          DO  ic = 1, pr_ncc
1871             SELECT CASE ( TRIM( cc_name(ic) ) )
1872                CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' )
1873                   cc_in2mod(1) = ic
1874                CASE ( 'OC', 'oc' )
1875                   cc_in2mod(2) = ic
1876                CASE ( 'BC', 'bc' )
1877                   cc_in2mod(3) = ic
1878                CASE ( 'DU', 'du' )
1879                   cc_in2mod(4) = ic
1880                CASE ( 'SS', 'ss' )
1881                   cc_in2mod(5) = ic
1882                CASE ( 'HNO3', 'hno3', 'NO3', 'no3', 'NO', 'no' )
1883                   cc_in2mod(6) = ic
1884                CASE ( 'NH3', 'nh3', 'NH4', 'nh4', 'NH', 'nh' )
1885                   cc_in2mod(7) = ic
1886             END SELECT
1887          ENDDO
1888
1889          IF ( SUM( cc_in2mod ) == 0 )  THEN
1890             message_string = 'None of the aerosol chemical components in ' // TRIM(               &
1891                              input_file_dynamic ) // ' correspond to ones applied in SALSA.'
1892             CALL message( 'salsa_mod: aerosol_init', 'PA0602', 2, 2, 0, 6, 0 )
1893          ENDIF
1894!
1895!--       Vertical profiles of mass fractions of different chemical components:
1896          IF ( check_existence( var_names, 'init_atmosphere_mass_fracs_a' ) )  THEN
1897             CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a,           &
1898                                0, pr_ncc-1, 0, pr_nz-1 )
1899          ELSE
1900             WRITE( message_string, * ) 'Missing init_atmosphere_mass_fracs_a in ' //              &
1901                                        TRIM( input_file_dynamic )
1902             CALL message( 'aerosol_init', 'PA0656', 1, 2, 0, 6, 0 )
1903          ENDIF
1904          CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_b', pr_mass_fracs_b,              &
1905                             0, pr_ncc-1, 0, pr_nz-1  )
1906!
1907!--       Match the input data with the chemical composition applied in the model
1908          DO  ic = 1, maxspec
1909             ss = cc_in2mod(ic)
1910             IF ( ss == 0 )  CYCLE
1911             pmf2a(nzb+1:nzt+1,ic) = pr_mass_fracs_a(nzb:nzt,ss)
1912             pmf2b(nzb+1:nzt+1,ic) = pr_mass_fracs_b(nzb:nzt,ss)
1913          ENDDO
1914!
1915!--       Aerosol concentrations: lod=1 (vertical profile of sectional number size distribution)
1916          CALL get_attribute( id_dyn, 'lod', lod_aero, .FALSE., 'init_atmosphere_aerosol' )
1917          IF ( lod_aero /= 1 )  THEN
1918             message_string = 'Currently only lod=1 accepted for init_atmosphere_aerosol'
1919             CALL message( 'salsa_mod: aerosol_init', 'PA0603', 2, 2, 0, 6, 0 )
1920          ELSE
1921!
1922!--          Bin mean diameters in the input file
1923             CALL get_dimension_length( id_dyn, pr_nbins, 'Dmid')
1924             IF ( pr_nbins /= nbins_aerosol )  THEN
1925                message_string = 'Number of size bins in init_atmosphere_aerosol does not match '  &
1926                                 // 'with that applied in the model'
1927                CALL message( 'salsa_mod: aerosol_init', 'PA0604', 2, 2, 0, 6, 0 )
1928             ENDIF
1929
1930             ALLOCATE( pr_dmid(pr_nbins) )
1931             pr_dmid    = 0.0_wp
1932
1933             CALL get_variable( id_dyn, 'Dmid', pr_dmid )
1934!
1935!--          Check whether the sectional representation conform to the one
1936!--          applied in the model
1937             IF ( ANY( ABS( ( aero(1:nbins_aerosol)%dmid - pr_dmid ) /                             &
1938                              aero(1:nbins_aerosol)%dmid )  > 0.1_wp )  ) THEN
1939                message_string = 'Mean diameters of the aerosol size bins in ' // TRIM(            &
1940                                 input_file_dynamic ) // ' do not match with the sectional '//     &
1941                                 'representation of the model.'
1942                CALL message( 'salsa_mod: aerosol_init', 'PA0605', 2, 2, 0, 6, 0 )
1943             ENDIF
1944!
1945!--          Inital aerosol concentrations
1946             CALL get_variable( id_dyn, 'init_atmosphere_aerosol', pndist(nzb+1:nzt,:),            &
1947                                0, pr_nbins-1, 0, pr_nz-1 )
1948          ENDIF
1949!
1950!--       Set bottom and top boundary condition (Neumann)
1951          pmf2a(nzb,:)    = pmf2a(nzb+1,:)
1952          pmf2a(nzt+1,:)  = pmf2a(nzt,:)
1953          pmf2b(nzb,:)    = pmf2b(nzb+1,:)
1954          pmf2b(nzt+1,:)  = pmf2b(nzt,:)
1955          pndist(nzb,:)   = pndist(nzb+1,:)
1956          pndist(nzt+1,:) = pndist(nzt,:)
1957
1958          IF ( index_so4 < 0 )  THEN
1959             pmf2a(:,1) = 0.0_wp
1960             pmf2b(:,1) = 0.0_wp
1961          ENDIF
1962          IF ( index_oc < 0 )  THEN
1963             pmf2a(:,2) = 0.0_wp
1964             pmf2b(:,2) = 0.0_wp
1965          ENDIF
1966          IF ( index_bc < 0 )  THEN
1967             pmf2a(:,3) = 0.0_wp
1968             pmf2b(:,3) = 0.0_wp
1969          ENDIF
1970          IF ( index_du < 0 )  THEN
1971             pmf2a(:,4) = 0.0_wp
1972             pmf2b(:,4) = 0.0_wp
1973          ENDIF
1974          IF ( index_ss < 0 )  THEN
1975             pmf2a(:,5) = 0.0_wp
1976             pmf2b(:,5) = 0.0_wp
1977          ENDIF
1978          IF ( index_no < 0 )  THEN
1979             pmf2a(:,6) = 0.0_wp
1980             pmf2b(:,6) = 0.0_wp
1981          ENDIF
1982          IF ( index_nh < 0 )  THEN
1983             pmf2a(:,7) = 0.0_wp
1984             pmf2b(:,7) = 0.0_wp
1985          ENDIF
1986
1987          IF ( SUM( pmf2a ) < 0.00001_wp  .AND.  SUM( pmf2b ) < 0.00001_wp )  THEN
1988             message_string = 'Error in initialising mass fractions of chemical components. ' //   &
1989                              'Check that all chemical components are included in parameter file!'
1990             CALL message( 'salsa_mod: aerosol_init', 'PA0606', 2, 2, 0, 6, 0 ) 
1991          ENDIF
1992!
1993!--       Then normalise the mass fraction so that SUM = 1
1994          DO  k = nzb, nzt+1
1995             pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) )
1996             IF ( SUM( pmf2b(k,:) ) > 0.0_wp )  pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) )
1997          ENDDO
1998
1999          DEALLOCATE( pr_z, pr_mass_fracs_a, pr_mass_fracs_b )
2000
2001       ELSE
2002          message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) //  &
2003                           ' for SALSA missing!'
2004          CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 )
2005!
2006!--       Close input file
2007          CALL close_input_file( id_dyn )
2008       ENDIF   ! netcdf_extend
2009
2010#else
2011       message_string = 'init_aerosol_type = 1 but preprocessor directive __netcdf is not used '// &
2012                        'in compiling!'
2013       CALL message( 'salsa_mod: aerosol_init', 'PA0608', 1, 2, 0, 6, 0 )
2014
2015#endif
2016
2017    ELSEIF ( init_aerosol_type == 0 )  THEN
2018!
2019!--    Mass fractions for species in a and b-bins
2020       IF ( index_so4 > 0 )  THEN
2021          pmf2a(:,1) = mass_fracs_a(index_so4)
2022          pmf2b(:,1) = mass_fracs_b(index_so4)
2023       ENDIF
2024       IF ( index_oc > 0 )  THEN
2025          pmf2a(:,2) = mass_fracs_a(index_oc)
2026          pmf2b(:,2) = mass_fracs_b(index_oc)
2027       ENDIF
2028       IF ( index_bc > 0 )  THEN
2029          pmf2a(:,3) = mass_fracs_a(index_bc)
2030          pmf2b(:,3) = mass_fracs_b(index_bc)
2031       ENDIF
2032       IF ( index_du > 0 )  THEN
2033          pmf2a(:,4) = mass_fracs_a(index_du)
2034          pmf2b(:,4) = mass_fracs_b(index_du)
2035       ENDIF
2036       IF ( index_ss > 0 )  THEN
2037          pmf2a(:,5) = mass_fracs_a(index_ss)
2038          pmf2b(:,5) = mass_fracs_b(index_ss)
2039       ENDIF
2040       IF ( index_no > 0 )  THEN
2041          pmf2a(:,6) = mass_fracs_a(index_no)
2042          pmf2b(:,6) = mass_fracs_b(index_no)
2043       ENDIF
2044       IF ( index_nh > 0 )  THEN
2045          pmf2a(:,7) = mass_fracs_a(index_nh)
2046          pmf2b(:,7) = mass_fracs_b(index_nh)
2047       ENDIF
2048       DO  k = nzb, nzt+1
2049          pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) )
2050          IF ( SUM( pmf2b(k,:) ) > 0.0_wp ) pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) )
2051       ENDDO
2052
2053       CALL size_distribution( n_lognorm, dpg, sigmag, nsect )
2054!
2055!--    Normalize by the given total number concentration
2056       nsect = nsect * SUM( n_lognorm ) / SUM( nsect )
2057       DO  ib = start_subrange_1a, end_subrange_2b
2058          pndist(:,ib) = nsect(ib)
2059       ENDDO
2060    ENDIF
2061
2062    IF ( init_gases_type == 1 )  THEN
2063!
2064!--    Read input profiles from PIDS_CHEM
2065#if defined( __netcdf )
2066!
2067!--    Location-dependent size distributions and compositions.
2068       INQUIRE( FILE = TRIM( input_file_dynamic ) //  TRIM( coupling_char ), EXIST = netcdf_extend )
2069       IF ( netcdf_extend  .AND.  .NOT. salsa_gases_from_chem )  THEN
2070!
2071!--       Open file in read-only mode
2072          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
2073!
2074!--       Inquire dimensions:
2075          CALL get_dimension_length( id_dyn, pr_nz, 'z' )
2076          IF ( pr_nz /= nz )  THEN
2077             WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//&
2078                                        'the number of numeric grid points.'
2079             CALL message( 'aerosol_init', 'PA0609', 1, 2, 0, 6, 0 )
2080          ENDIF
2081!
2082!--       Read vertical profiles of gases:
2083          CALL get_variable( id_dyn, 'init_atmosphere_h2so4', salsa_gas(1)%init(nzb+1:nzt) )
2084          CALL get_variable( id_dyn, 'init_atmosphere_hno3',  salsa_gas(2)%init(nzb+1:nzt) )
2085          CALL get_variable( id_dyn, 'init_atmosphere_nh3',   salsa_gas(3)%init(nzb+1:nzt) )
2086          CALL get_variable( id_dyn, 'init_atmosphere_ocnv',  salsa_gas(4)%init(nzb+1:nzt) )
2087          CALL get_variable( id_dyn, 'init_atmosphere_ocsv',  salsa_gas(5)%init(nzb+1:nzt) )
2088!
2089!--       Set Neumann top and surface boundary condition for initial + initialise concentrations
2090          DO  ig = 1, ngases_salsa
2091             salsa_gas(ig)%init(nzb)   =  salsa_gas(ig)%init(nzb+1)
2092             salsa_gas(ig)%init(nzt+1) =  salsa_gas(ig)%init(nzt)
2093             IF ( .NOT. read_restart_data_salsa )  THEN
2094                DO  k = nzb, nzt+1
2095                   salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
2096                ENDDO
2097             ENDIF
2098          ENDDO
2099
2100       ELSEIF ( .NOT. netcdf_extend  .AND.  .NOT.  salsa_gases_from_chem )  THEN
2101          message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) //  &
2102                           ' for SALSA missing!'
2103          CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 )
2104!
2105!--       Close input file
2106          CALL close_input_file( id_dyn )
2107       ENDIF   ! netcdf_extend
2108#else
2109       message_string = 'init_gases_type = 1 but preprocessor directive __netcdf is not used in '//&
2110                        'compiling!'
2111       CALL message( 'salsa_mod: aerosol_init', 'PA0611', 1, 2, 0, 6, 0 )
2112
2113#endif
2114
2115    ENDIF
2116!
2117!-- Both SO4 and OC are included, so use the given mass fractions
2118    IF ( index_oc > 0  .AND.  index_so4 > 0 )  THEN
2119       pmfoc1a(:) = pmf2a(:,2) / ( pmf2a(:,2) + pmf2a(:,1) )  ! Normalize
2120!
2121!-- Pure organic carbon
2122    ELSEIF ( index_oc > 0 )  THEN
2123       pmfoc1a(:) = 1.0_wp
2124!
2125!-- Pure SO4
2126    ELSEIF ( index_so4 > 0 )  THEN
2127       pmfoc1a(:) = 0.0_wp
2128
2129    ELSE
2130       message_string = 'Either OC or SO4 must be active for aerosol region 1a!'
2131       CALL message( 'salsa_mod: aerosol_init', 'PA0612', 1, 2, 0, 6, 0 )
2132    ENDIF
2133
2134!
2135!-- Initialize concentrations
2136    DO  i = nxlg, nxrg
2137       DO  j = nysg, nyng
2138          DO  k = nzb, nzt+1
2139!
2140!--          Predetermine flag to mask topography
2141             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2142!
2143!--          a) Number concentrations
2144!--          Region 1:
2145             DO  ib = start_subrange_1a, end_subrange_1a
2146                IF ( .NOT. read_restart_data_salsa )  THEN
2147                   aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
2148                ENDIF
2149                IF ( prunmode == 1 )  THEN
2150                   aerosol_number(ib)%init = pndist(:,ib)
2151                ENDIF
2152             ENDDO
2153!
2154!--          Region 2:
2155             IF ( nreg > 1 )  THEN
2156                DO  ib = start_subrange_2a, end_subrange_2a
2157                   IF ( .NOT. read_restart_data_salsa )  THEN
2158                      aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
2159                   ENDIF
2160                   IF ( prunmode == 1 )  THEN
2161                      aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib)
2162                   ENDIF
2163                ENDDO
2164                IF ( .NOT. no_insoluble )  THEN
2165                   DO  ib = start_subrange_2b, end_subrange_2b
2166                      IF ( pnf2a(k) < 1.0_wp )  THEN
2167                         IF ( .NOT. read_restart_data_salsa )  THEN
2168                            aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *    &
2169                                                             pndist(k,ib) * flag
2170                         ENDIF
2171                         IF ( prunmode == 1 )  THEN
2172                            aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp - nf2a ) * pndist(:,ib)
2173                         ENDIF
2174                      ENDIF
2175                   ENDDO
2176                ENDIF
2177             ENDIF
2178!
2179!--          b) Aerosol mass concentrations
2180!--             bin subrange 1: done here separately due to the SO4/OC convention
2181!
2182!--          SO4:
2183             IF ( index_so4 > 0 )  THEN
2184                ss = ( index_so4 - 1 ) * nbins_aerosol + start_subrange_1a !< start
2185                ee = ( index_so4 - 1 ) * nbins_aerosol + end_subrange_1a !< end
2186                ib = start_subrange_1a
2187                DO  ic = ss, ee
2188                   IF ( .NOT. read_restart_data_salsa )  THEN
2189                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) *          &
2190                                                     pndist(k,ib) * core(ib) * arhoh2so4 * flag
2191                   ENDIF
2192                   IF ( prunmode == 1 )  THEN
2193                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib) &
2194                                                 * core(ib) * arhoh2so4
2195                   ENDIF
2196                   ib = ib+1
2197                ENDDO
2198             ENDIF
2199!
2200!--          OC:
2201             IF ( index_oc > 0 ) THEN
2202                ss = ( index_oc - 1 ) * nbins_aerosol + start_subrange_1a !< start
2203                ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a !< end
2204                ib = start_subrange_1a
2205                DO  ic = ss, ee
2206                   IF ( .NOT. read_restart_data_salsa )  THEN
2207                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *    &
2208                                                     core(ib) * arhooc * flag
2209                   ENDIF
2210                   IF ( prunmode == 1 )  THEN
2211                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *        &
2212                                                 core(ib) * arhooc
2213                   ENDIF
2214                   ib = ib+1
2215                ENDDO
2216             ENDIF
2217          ENDDO !< k
2218
2219          prunmode = 3  ! Init only once
2220
2221       ENDDO !< j
2222    ENDDO !< i
2223
2224!
2225!-- c) Aerosol mass concentrations
2226!--    bin subrange 2:
2227    IF ( nreg > 1 ) THEN
2228
2229       IF ( index_so4 > 0 ) THEN
2230          CALL set_aero_mass( index_so4, pmf2a(:,1), pmf2b(:,1), pnf2a, pndist, core, arhoh2so4 )
2231       ENDIF
2232       IF ( index_oc > 0 ) THEN
2233          CALL set_aero_mass( index_oc, pmf2a(:,2), pmf2b(:,2), pnf2a, pndist, core, arhooc )
2234       ENDIF
2235       IF ( index_bc > 0 ) THEN
2236          CALL set_aero_mass( index_bc, pmf2a(:,3), pmf2b(:,3), pnf2a, pndist, core, arhobc )
2237       ENDIF
2238       IF ( index_du > 0 ) THEN
2239          CALL set_aero_mass( index_du, pmf2a(:,4), pmf2b(:,4), pnf2a, pndist, core, arhodu )
2240       ENDIF
2241       IF ( index_ss > 0 ) THEN
2242          CALL set_aero_mass( index_ss, pmf2a(:,5), pmf2b(:,5), pnf2a, pndist, core, arhoss )
2243       ENDIF
2244       IF ( index_no > 0 ) THEN
2245          CALL set_aero_mass( index_no, pmf2a(:,6), pmf2b(:,6), pnf2a, pndist, core, arhohno3 )
2246       ENDIF
2247       IF ( index_nh > 0 ) THEN
2248          CALL set_aero_mass( index_nh, pmf2a(:,7), pmf2b(:,7), pnf2a, pndist, core, arhonh3 )
2249       ENDIF
2250
2251    ENDIF
2252
2253 END SUBROUTINE aerosol_init
2254
2255!------------------------------------------------------------------------------!
2256! Description:
2257! ------------
2258!> Create a lognormal size distribution and discretise to a sectional
2259!> representation.
2260!------------------------------------------------------------------------------!
2261 SUBROUTINE size_distribution( in_ntot, in_dpg, in_sigma, psd_sect )
2262
2263    IMPLICIT NONE
2264
2265    INTEGER(iwp) ::  ib         !< running index: bin
2266    INTEGER(iwp) ::  iteration  !< running index: iteration
2267
2268    REAL(wp) ::  d1         !< particle diameter (m, dummy)
2269    REAL(wp) ::  d2         !< particle diameter (m, dummy)
2270    REAL(wp) ::  delta_d    !< (d2-d1)/10
2271    REAL(wp) ::  deltadp    !< bin width
2272    REAL(wp) ::  dmidi      !< ( d1 + d2 ) / 2
2273
2274    REAL(wp), DIMENSION(:), INTENT(in) ::  in_dpg    !< geometric mean diameter (m)
2275    REAL(wp), DIMENSION(:), INTENT(in) ::  in_ntot   !< number conc. (#/m3)
2276    REAL(wp), DIMENSION(:), INTENT(in) ::  in_sigma  !< standard deviation
2277
2278    REAL(wp), DIMENSION(:), INTENT(inout) ::  psd_sect  !< sectional size distribution
2279
2280    DO  ib = start_subrange_1a, end_subrange_2b
2281       psd_sect(ib) = 0.0_wp
2282!
2283!--    Particle diameter at the low limit (largest in the bin) (m)
2284       d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp
2285!
2286!--    Particle diameter at the high limit (smallest in the bin) (m)
2287       d2 = ( aero(ib)%vhilim / api6 )**0.33333333_wp
2288!
2289!--    Span of particle diameter in a bin (m)
2290       delta_d = 0.1_wp * ( d2 - d1 )
2291!
2292!--    Iterate:
2293       DO  iteration = 1, 10
2294          d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp + ( ib - 1) * delta_d
2295          d2 = d1 + delta_d
2296          dmidi = 0.5_wp * ( d1 + d2 )
2297          deltadp = LOG10( d2 / d1 )
2298!
2299!--       Size distribution
2300!--       in_ntot = total number, total area, or total volume concentration
2301!--       in_dpg = geometric-mean number, area, or volume diameter
2302!--       n(k) = number, area, or volume concentration in a bin
2303          psd_sect(ib) = psd_sect(ib) + SUM( in_ntot * deltadp / ( SQRT( 2.0_wp * pi ) *           &
2304                        LOG10( in_sigma ) ) * EXP( -LOG10( dmidi / in_dpg )**2.0_wp /              &
2305                        ( 2.0_wp * LOG10( in_sigma ) ** 2.0_wp ) ) )
2306
2307       ENDDO
2308    ENDDO
2309
2310 END SUBROUTINE size_distribution
2311
2312!------------------------------------------------------------------------------!
2313! Description:
2314! ------------
2315!> Sets the mass concentrations to aerosol arrays in 2a and 2b.
2316!>
2317!> Tomi Raatikainen, FMI, 29.2.2016
2318!------------------------------------------------------------------------------!
2319 SUBROUTINE set_aero_mass( ispec, pmf2a, pmf2b, pnf2a, pndist, pcore, prho )
2320
2321    IMPLICIT NONE
2322
2323    INTEGER(iwp) ::  ee        !< index: end
2324    INTEGER(iwp) ::  i         !< loop index
2325    INTEGER(iwp) ::  ib        !< loop index
2326    INTEGER(iwp) ::  ic        !< loop index
2327    INTEGER(iwp) ::  j         !< loop index
2328    INTEGER(iwp) ::  k         !< loop index
2329    INTEGER(iwp) ::  prunmode  !< 1 = initialise
2330    INTEGER(iwp) ::  ss        !< index: start
2331
2332    INTEGER(iwp), INTENT(in) :: ispec  !< Aerosol species index
2333
2334    REAL(wp) ::  flag   !< flag to mask topography grid points
2335
2336    REAL(wp), INTENT(in) ::  prho !< Aerosol density
2337
2338    REAL(wp), DIMENSION(nbins_aerosol), INTENT(in) ::  pcore !< Aerosol bin mid core volume
2339    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pnf2a !< Number fraction for 2a
2340    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pmf2a !< Mass distributions for a
2341    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pmf2b !< and b bins
2342
2343    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol), INTENT(in) ::  pndist !< Aerosol size distribution
2344
2345    prunmode = 1
2346
2347    DO i = nxlg, nxrg
2348       DO j = nysg, nyng
2349          DO k = nzb, nzt+1
2350!
2351!--          Predetermine flag to mask topography
2352             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
2353!
2354!--          Regime 2a:
2355             ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2a
2356             ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2a
2357             ib = start_subrange_2a
2358             DO ic = ss, ee
2359                IF ( .NOT. read_restart_data_salsa )  THEN
2360                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib)&
2361                                                  * pcore(ib) * prho * flag
2362                ENDIF
2363                IF ( prunmode == 1 )  THEN
2364                   aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) *  &
2365                                              pcore(ib) * prho
2366                ENDIF
2367                ib = ib + 1
2368             ENDDO
2369!
2370!--          Regime 2b:
2371             IF ( .NOT. no_insoluble )  THEN
2372                ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2b
2373                ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2b
2374                ib = start_subrange_2a
2375                DO ic = ss, ee
2376                   IF ( .NOT. read_restart_data_salsa )  THEN
2377                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k))&
2378                                                     * pndist(k,ib) * pcore(ib) * prho * flag
2379                   ENDIF
2380                   IF ( prunmode == 1 )  THEN
2381                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) * &
2382                                                 pndist(k,ib) * pcore(ib) * prho
2383                   ENDIF
2384                   ib = ib + 1
2385                ENDDO  ! c
2386
2387             ENDIF
2388          ENDDO   ! k
2389
2390          prunmode = 3  ! Init only once
2391
2392       ENDDO   ! j
2393    ENDDO   ! i
2394
2395 END SUBROUTINE set_aero_mass
2396
2397!------------------------------------------------------------------------------!
2398! Description:
2399! ------------
2400!> Initialise the matching between surface types in LSM and deposition models.
2401!> Do the matching based on Zhang et al. (2001). Atmos. Environ. 35, 549-560
2402!> (here referred as Z01).
2403!------------------------------------------------------------------------------!
2404 SUBROUTINE init_deposition
2405
2406    USE surface_mod,                                                                               &
2407        ONLY:  surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
2408
2409    IMPLICIT NONE
2410
2411    INTEGER(iwp) ::  l  !< loop index for vertical surfaces
2412
2413    LOGICAL :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM surfaces)
2414
2415    IF ( depo_pcm_par == 'zhang2001' )  THEN
2416       depo_pcm_par_num = 1
2417    ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
2418       depo_pcm_par_num = 2
2419    ENDIF
2420
2421    IF ( depo_surf_par == 'zhang2001' )  THEN
2422       depo_surf_par_num = 1
2423    ELSEIF ( depo_surf_par == 'petroff2010' )  THEN
2424       depo_surf_par_num = 2
2425    ENDIF
2426!
2427!-- LSM: Pavement, vegetation and water
2428    IF ( nldepo_surf  .AND.  land_surface )  THEN
2429       match_lsm = .TRUE.
2430       ALLOCATE( lsm_to_depo_h%match_lupg(1:surf_lsm_h%ns),                                         &
2431                 lsm_to_depo_h%match_luvw(1:surf_lsm_h%ns),                                         &
2432                 lsm_to_depo_h%match_luww(1:surf_lsm_h%ns) )
2433       lsm_to_depo_h%match_lupg = 0
2434       lsm_to_depo_h%match_luvw = 0
2435       lsm_to_depo_h%match_luww = 0
2436       CALL match_sm_zhang( surf_lsm_h, lsm_to_depo_h%match_lupg, lsm_to_depo_h%match_luvw,        &
2437                            lsm_to_depo_h%match_luww, match_lsm )
2438       DO  l = 0, 3
2439          ALLOCATE( lsm_to_depo_v(l)%match_lupg(1:surf_lsm_v(l)%ns),                               &
2440                    lsm_to_depo_v(l)%match_luvw(1:surf_lsm_v(l)%ns),                               &
2441                    lsm_to_depo_v(l)%match_luww(1:surf_lsm_v(l)%ns) )
2442          lsm_to_depo_v(l)%match_lupg = 0
2443          lsm_to_depo_v(l)%match_luvw = 0
2444          lsm_to_depo_v(l)%match_luww = 0
2445          CALL match_sm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match_lupg,                         &
2446                               lsm_to_depo_v(l)%match_luvw, lsm_to_depo_v(l)%match_luww, match_lsm )
2447       ENDDO
2448    ENDIF
2449!
2450!-- USM: Green roofs/walls, wall surfaces and windows
2451    IF ( nldepo_surf  .AND.  urban_surface )  THEN
2452       match_lsm = .FALSE.
2453       ALLOCATE( usm_to_depo_h%match_lupg(1:surf_usm_h%ns),                                        &
2454                 usm_to_depo_h%match_luvw(1:surf_usm_h%ns),                                        &
2455                 usm_to_depo_h%match_luww(1:surf_usm_h%ns) )
2456       usm_to_depo_h%match_lupg = 0
2457       usm_to_depo_h%match_luvw = 0
2458       usm_to_depo_h%match_luww = 0
2459       CALL match_sm_zhang( surf_usm_h, usm_to_depo_h%match_lupg, usm_to_depo_h%match_luvw,        &
2460                            usm_to_depo_h%match_luww, match_lsm )
2461       DO  l = 0, 3
2462          ALLOCATE( usm_to_depo_v(l)%match_lupg(1:surf_usm_v(l)%ns),                               &
2463                    usm_to_depo_v(l)%match_luvw(1:surf_usm_v(l)%ns),                               &
2464                    usm_to_depo_v(l)%match_luww(1:surf_usm_v(l)%ns) )
2465          usm_to_depo_v(l)%match_lupg = 0
2466          usm_to_depo_v(l)%match_luvw = 0
2467          usm_to_depo_v(l)%match_luww = 0
2468          CALL match_sm_zhang( surf_usm_v(l), usm_to_depo_v(l)%match_lupg,                         &
2469                               usm_to_depo_v(l)%match_luvw, usm_to_depo_v(l)%match_luww, match_lsm )
2470       ENDDO
2471    ENDIF
2472
2473    IF ( nldepo_pcm )  THEN
2474       SELECT CASE ( depo_pcm_type )
2475          CASE ( 'evergreen_needleleaf' )
2476             depo_pcm_type_num = 1
2477          CASE ( 'evergreen_broadleaf' )
2478             depo_pcm_type_num = 2
2479          CASE ( 'deciduous_needleleaf' )
2480             depo_pcm_type_num = 3
2481          CASE ( 'deciduous_broadleaf' )
2482             depo_pcm_type_num = 4
2483          CASE DEFAULT
2484             message_string = 'depo_pcm_type not set correctly.'
2485             CALL message( 'salsa_mod: init_deposition', 'PA0613', 1, 2, 0, 6, 0 )
2486       END SELECT
2487    ENDIF
2488
2489 END SUBROUTINE init_deposition
2490
2491!------------------------------------------------------------------------------!
2492! Description:
2493! ------------
2494!> Match the surface types in PALM and Zhang et al. 2001 deposition module
2495!------------------------------------------------------------------------------!
2496 SUBROUTINE match_sm_zhang( surf, match_pav_green, match_veg_wall, match_wat_win, match_lsm )
2497
2498    USE surface_mod,                                                           &
2499        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
2500
2501    IMPLICIT NONE
2502
2503    INTEGER(iwp) ::  m              !< index for surface elements
2504    INTEGER(iwp) ::  pav_type_palm  !< pavement / green wall type in PALM
2505    INTEGER(iwp) ::  veg_type_palm  !< vegetation / wall type in PALM
2506    INTEGER(iwp) ::  wat_type_palm  !< water / window type in PALM
2507
2508    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_pav_green  !<  matching pavement/green walls
2509    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_veg_wall   !<  matching vegetation/walls
2510    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_wat_win    !<  matching water/windows
2511
2512    LOGICAL, INTENT(in) :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM)
2513
2514    TYPE(surf_type), INTENT(in) :: surf  !< respective surface type
2515
2516    DO  m = 1, surf%ns
2517       IF ( match_lsm )  THEN
2518!
2519!--       Vegetation (LSM):
2520          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
2521             veg_type_palm = surf%vegetation_type(m)
2522             SELECT CASE ( veg_type_palm )
2523                CASE ( 0 )
2524                   message_string = 'No vegetation type defined.'
2525                   CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
2526                CASE ( 1 )  ! bare soil
2527                   match_veg_wall(m) = 6  ! grass in Z01
2528                CASE ( 2 )  ! crops, mixed farming
2529                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
2530                CASE ( 3 )  ! short grass
2531                   match_veg_wall(m) = 6  ! grass in Z01
2532                CASE ( 4 )  ! evergreen needleleaf trees
2533                    match_veg_wall(m) = 1  ! evergreen needleleaf trees in Z01
2534                CASE ( 5 )  ! deciduous needleleaf trees
2535                   match_veg_wall(m) = 3  ! deciduous needleleaf trees in Z01
2536                CASE ( 6 )  ! evergreen broadleaf trees
2537                   match_veg_wall(m) = 2  ! evergreen broadleaf trees in Z01
2538                CASE ( 7 )  ! deciduous broadleaf trees
2539                   match_veg_wall(m) = 4  ! deciduous broadleaf trees in Z01
2540                CASE ( 8 )  ! tall grass
2541                   match_veg_wall(m) = 6  ! grass in Z01
2542                CASE ( 9 )  ! desert
2543                   match_veg_wall(m) = 8  ! desert in Z01
2544                CASE ( 10 )  ! tundra
2545                   match_veg_wall(m) = 9  ! tundra in Z01
2546                CASE ( 11 )  ! irrigated crops
2547                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
2548                CASE ( 12 )  ! semidesert
2549                   match_veg_wall(m) = 8  ! desert in Z01
2550                CASE ( 13 )  ! ice caps and glaciers
2551                   match_veg_wall(m) = 12  ! ice cap and glacier in Z01
2552                CASE ( 14 )  ! bogs and marshes
2553                   match_veg_wall(m) = 11  ! wetland with plants in Z01
2554                CASE ( 15 )  ! evergreen shrubs
2555                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2556                CASE ( 16 )  ! deciduous shrubs
2557                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2558                CASE ( 17 )  ! mixed forest/woodland
2559                   match_veg_wall(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
2560                CASE ( 18 )  ! interrupted forest
2561                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2562             END SELECT
2563          ENDIF
2564!
2565!--       Pavement (LSM):
2566          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
2567             pav_type_palm = surf%pavement_type(m)
2568             IF ( pav_type_palm == 0 )  THEN  ! error
2569                message_string = 'No pavement type defined.'
2570                CALL message( 'salsa_mod: match_sm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
2571             ELSE
2572                match_pav_green(m) = 15  ! urban in Z01
2573             ENDIF
2574          ENDIF
2575!
2576!--       Water (LSM):
2577          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
2578             wat_type_palm = surf%water_type(m)
2579             IF ( wat_type_palm == 0 )  THEN  ! error
2580                message_string = 'No water type defined.'
2581                CALL message( 'salsa_mod: match_sm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
2582             ELSEIF ( wat_type_palm == 3 )  THEN
2583                match_wat_win(m) = 14  ! ocean in Z01
2584             ELSEIF ( wat_type_palm == 1  .OR.  wat_type_palm == 2 .OR.  wat_type_palm == 4        &
2585                      .OR.  wat_type_palm == 5  )  THEN
2586                match_wat_win(m) = 13  ! inland water in Z01
2587             ENDIF
2588          ENDIF
2589       ELSE
2590!
2591!--       Wall surfaces (USM):
2592          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
2593             match_veg_wall(m) = 15  ! urban in Z01
2594          ENDIF
2595!
2596!--       Green walls and roofs (USM):
2597          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
2598             match_pav_green(m) =  6 ! (short) grass in Z01
2599          ENDIF
2600!
2601!--       Windows (USM):
2602          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
2603             match_wat_win(m) = 15  ! urban in Z01
2604          ENDIF
2605       ENDIF
2606
2607    ENDDO
2608
2609 END SUBROUTINE match_sm_zhang
2610
2611!------------------------------------------------------------------------------!
2612! Description:
2613! ------------
2614!> Swapping of timelevels
2615!------------------------------------------------------------------------------!
2616 SUBROUTINE salsa_swap_timelevel( mod_count )
2617
2618    IMPLICIT NONE
2619
2620    INTEGER(iwp) ::  ib   !<
2621    INTEGER(iwp) ::  ic   !<
2622    INTEGER(iwp) ::  icc  !<
2623    INTEGER(iwp) ::  ig   !<
2624
2625    INTEGER(iwp), INTENT(IN) ::  mod_count  !<
2626
2627    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
2628
2629       SELECT CASE ( mod_count )
2630
2631          CASE ( 0 )
2632
2633             DO  ib = 1, nbins_aerosol
2634                aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => nconc_1(:,:,:,ib)
2635                aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,ib)
2636
2637                DO  ic = 1, ncomponents_mass
2638                   icc = ( ic-1 ) * nbins_aerosol + ib
2639                   aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => mconc_1(:,:,:,icc)
2640                   aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,icc)
2641                ENDDO
2642             ENDDO
2643
2644             IF ( .NOT. salsa_gases_from_chem )  THEN
2645                DO  ig = 1, ngases_salsa
2646                   salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => gconc_1(:,:,:,ig)
2647                   salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,ig)
2648                ENDDO
2649             ENDIF
2650
2651          CASE ( 1 )
2652
2653             DO  ib = 1, nbins_aerosol
2654                aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => nconc_2(:,:,:,ib)
2655                aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,ib)
2656                DO  ic = 1, ncomponents_mass
2657                   icc = ( ic-1 ) * nbins_aerosol + ib
2658                   aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => mconc_2(:,:,:,icc)
2659                   aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,icc)
2660                ENDDO
2661             ENDDO
2662
2663             IF ( .NOT. salsa_gases_from_chem )  THEN
2664                DO  ig = 1, ngases_salsa
2665                   salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => gconc_2(:,:,:,ig)
2666                   salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,ig)
2667                ENDDO
2668             ENDIF
2669
2670       END SELECT
2671
2672    ENDIF
2673
2674 END SUBROUTINE salsa_swap_timelevel
2675
2676
2677!------------------------------------------------------------------------------!
2678! Description:
2679! ------------
2680!> This routine reads the respective restart data.
2681!------------------------------------------------------------------------------!
2682 SUBROUTINE salsa_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,      &
2683                             nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
2684
2685    USE control_parameters,                                                                        &
2686        ONLY:  length, restart_string
2687
2688    IMPLICIT NONE
2689
2690    INTEGER(iwp) ::  ib              !<
2691    INTEGER(iwp) ::  ic              !<
2692    INTEGER(iwp) ::  ig              !<
2693    INTEGER(iwp) ::  k               !<
2694    INTEGER(iwp) ::  nxlc            !<
2695    INTEGER(iwp) ::  nxlf            !<
2696    INTEGER(iwp) ::  nxl_on_file     !<
2697    INTEGER(iwp) ::  nxrc            !<
2698    INTEGER(iwp) ::  nxrf            !<
2699    INTEGER(iwp) ::  nxr_on_file     !<
2700    INTEGER(iwp) ::  nync            !<
2701    INTEGER(iwp) ::  nynf            !<
2702    INTEGER(iwp) ::  nyn_on_file     !<
2703    INTEGER(iwp) ::  nysc            !<
2704    INTEGER(iwp) ::  nysf            !<
2705    INTEGER(iwp) ::  nys_on_file     !<
2706
2707    LOGICAL, INTENT(OUT)  ::  found  !<
2708
2709    REAL(wp), &
2710       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
2711
2712    found = .FALSE.
2713
2714    IF ( read_restart_data_salsa )  THEN
2715
2716       SELECT CASE ( restart_string(1:length) )
2717
2718          CASE ( 'aerosol_number' )
2719             DO  ib = 1, nbins_aerosol
2720                IF ( k == 1 )  READ ( 13 ) tmp_3d
2721                aerosol_number(ib)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
2722                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2723                found = .TRUE.
2724             ENDDO
2725
2726          CASE ( 'aerosol_mass' )
2727             DO  ic = 1, ncomponents_mass * nbins_aerosol
2728                IF ( k == 1 )  READ ( 13 ) tmp_3d
2729                aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
2730                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2731                found = .TRUE.
2732             ENDDO
2733
2734          CASE ( 'salsa_gas' )
2735             DO  ig = 1, ngases_salsa
2736                IF ( k == 1 )  READ ( 13 ) tmp_3d
2737                salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                    &
2738                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2739                found = .TRUE.
2740             ENDDO
2741
2742          CASE DEFAULT
2743             found = .FALSE.
2744
2745       END SELECT
2746    ENDIF
2747
2748 END SUBROUTINE salsa_rrd_local
2749
2750!------------------------------------------------------------------------------!
2751! Description:
2752! ------------
2753!> This routine writes the respective restart data.
2754!> Note that the following input variables in PARIN have to be equal between
2755!> restart runs:
2756!>    listspec, nbin, nbin2, nf2a, ncc, mass_fracs_a, mass_fracs_b
2757!------------------------------------------------------------------------------!
2758 SUBROUTINE salsa_wrd_local
2759
2760    USE control_parameters,                                                                        &
2761        ONLY:  write_binary
2762
2763    IMPLICIT NONE
2764
2765    INTEGER(iwp) ::  ib   !<
2766    INTEGER(iwp) ::  ic   !<
2767    INTEGER(iwp) ::  ig  !<
2768
2769    IF ( write_binary  .AND.  write_binary_salsa )  THEN
2770
2771       CALL wrd_write_string( 'aerosol_number' )
2772       DO  ib = 1, nbins_aerosol
2773          WRITE ( 14 )  aerosol_number(ib)%conc
2774       ENDDO
2775
2776       CALL wrd_write_string( 'aerosol_mass' )
2777       DO  ic = 1, nbins_aerosol * ncomponents_mass
2778          WRITE ( 14 )  aerosol_mass(ic)%conc
2779       ENDDO
2780
2781       CALL wrd_write_string( 'salsa_gas' )
2782       DO  ig = 1, ngases_salsa
2783          WRITE ( 14 )  salsa_gas(ig)%conc
2784       ENDDO
2785
2786    ENDIF
2787
2788 END SUBROUTINE salsa_wrd_local
2789
2790!------------------------------------------------------------------------------!
2791! Description:
2792! ------------
2793!> Performs necessary unit and dimension conversion between the host model and
2794!> SALSA module, and calls the main SALSA routine.
2795!> Partially adobted form the original SALSA boxmodel version.
2796!> Now takes masses in as kg/kg from LES!! Converted to m3/m3 for SALSA
2797!> 05/2016 Juha: This routine is still pretty much in its original shape.
2798!>               It's dumb as a mule and twice as ugly, so implementation of
2799!>               an improved solution is necessary sooner or later.
2800!> Juha Tonttila, FMI, 2014
2801!> Jaakko Ahola, FMI, 2016
2802!> Only aerosol processes included, Mona Kurppa, UHel, 2017
2803!------------------------------------------------------------------------------!
2804 SUBROUTINE salsa_driver( i, j, prunmode )
2805
2806    USE arrays_3d,                                                                                 &
2807        ONLY: pt_p, q_p, u, v, w
2808
2809    USE plant_canopy_model_mod,                                                                    &
2810        ONLY: lad_s
2811
2812    USE surface_mod,                                                                               &
2813        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
2814
2815    IMPLICIT NONE
2816
2817    INTEGER(iwp) ::  endi    !< end index
2818    INTEGER(iwp) ::  ib      !< loop index
2819    INTEGER(iwp) ::  ic      !< loop index
2820    INTEGER(iwp) ::  ig      !< loop index
2821    INTEGER(iwp) ::  k_wall  !< vertical index of topography top
2822    INTEGER(iwp) ::  k       !< loop index
2823    INTEGER(iwp) ::  l       !< loop index
2824    INTEGER(iwp) ::  nc_h2o  !< index of H2O in the prtcl index table
2825    INTEGER(iwp) ::  ss      !< loop index
2826    INTEGER(iwp) ::  str     !< start index
2827    INTEGER(iwp) ::  vc      !< default index in prtcl
2828
2829    INTEGER(iwp), INTENT(in) ::  i         !< loop index
2830    INTEGER(iwp), INTENT(in) ::  j         !< loop index
2831    INTEGER(iwp), INTENT(in) ::  prunmode  !< 1: Initialization, 2: Spinup, 3: Regular runtime
2832
2833    REAL(wp) ::  cw_old  !< previous H2O mixing ratio
2834    REAL(wp) ::  flag    !< flag to mask topography grid points
2835    REAL(wp) ::  in_lad  !< leaf area density (m2/m3)
2836    REAL(wp) ::  in_rh   !< relative humidity
2837    REAL(wp) ::  zgso4   !< SO4
2838    REAL(wp) ::  zghno3  !< HNO3
2839    REAL(wp) ::  zgnh3   !< NH3
2840    REAL(wp) ::  zgocnv  !< non-volatile OC
2841    REAL(wp) ::  zgocsv  !< semi-volatile OC
2842
2843    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_adn  !< air density (kg/m3)
2844    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cs   !< H2O sat. vapour conc.
2845    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cw   !< H2O vapour concentration
2846    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_p    !< pressure (Pa)
2847    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_t    !< temperature (K)
2848    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_u    !< wind magnitude (m/s)
2849    REAL(wp), DIMENSION(nzb:nzt+1) ::  kvis    !< kinematic viscosity of air(m2/s)
2850    REAL(wp), DIMENSION(nzb:nzt+1) ::  ppm_to_nconc  !< Conversion factor from ppm to #/m3
2851
2852    REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) ::  schmidt_num  !< particle Schmidt number
2853    REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) ::  vd           !< particle fall seed (m/s)
2854
2855    TYPE(t_section), DIMENSION(nbins_aerosol) ::  lo_aero   !< additional variable for OpenMP
2856    TYPE(t_section), DIMENSION(nbins_aerosol) ::  aero_old  !< helper array
2857
2858    aero_old(:)%numc = 0.0_wp
2859    in_lad           = 0.0_wp
2860    in_u             = 0.0_wp
2861    kvis             = 0.0_wp
2862    lo_aero          = aero
2863    schmidt_num      = 0.0_wp
2864    vd               = 0.0_wp
2865    zgso4            = nclim
2866    zghno3           = nclim
2867    zgnh3            = nclim
2868    zgocnv           = nclim
2869    zgocsv           = nclim
2870!
2871!-- Aerosol number is always set, but mass can be uninitialized
2872    DO ib = 1, nbins_aerosol
2873       lo_aero(ib)%volc(:)  = 0.0_wp
2874       aero_old(ib)%volc(:) = 0.0_wp
2875    ENDDO
2876!
2877!-- Set the salsa runtime config (How to make this more efficient?)
2878    CALL set_salsa_runtime( prunmode )
2879!
2880!-- Calculate thermodynamic quantities needed in SALSA
2881    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 )
2882!
2883!-- Magnitude of wind: needed for deposition
2884    IF ( lsdepo )  THEN
2885       in_u(nzb+1:nzt) = SQRT( ( 0.5_wp * ( u(nzb+1:nzt,j,i) + u(nzb+1:nzt,j,i+1) ) )**2 +         &
2886                               ( 0.5_wp * ( v(nzb+1:nzt,j,i) + v(nzb+1:nzt,j+1,i) ) )**2 +         &
2887                               ( 0.5_wp * ( w(nzb:nzt-1,j,i) + w(nzb+1:nzt,j,  i) ) )**2 )
2888    ENDIF
2889!
2890!-- Calculate conversion factors for gas concentrations
2891    ppm_to_nconc(:) = for_ppm_to_nconc * in_p(:) / in_t(:)
2892!
2893!-- Determine topography-top index on scalar grid
2894    k_wall = k_topo_top(j,i)
2895
2896    DO k = nzb+1, nzt
2897!
2898!--    Predetermine flag to mask topography
2899       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2900!
2901!--    Wind velocity for dry depositon on vegetation
2902       IF ( lsdepo_pcm  .AND.  plant_canopy )  THEN
2903          in_lad = lad_s( MAX( k-k_wall,0 ),j,i)
2904       ENDIF
2905!
2906!--    For initialization and spinup, limit the RH with the parameter rhlim
2907       IF ( prunmode < 3 ) THEN
2908          in_cw(k) = MIN( in_cw(k), in_cs(k) * rhlim )
2909       ELSE
2910          in_cw(k) = in_cw(k)
2911       ENDIF
2912       cw_old = in_cw(k) !* in_adn(k)
2913!
2914!--    Set volume concentrations:
2915!--    Sulphate (SO4) or sulphuric acid H2SO4
2916       IF ( index_so4 > 0 )  THEN
2917          vc = 1
2918          str = ( index_so4-1 ) * nbins_aerosol + 1    ! start index
2919          endi = index_so4 * nbins_aerosol             ! end index
2920          ic = 1
2921          DO ss = str, endi
2922             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2so4
2923             ic = ic+1
2924          ENDDO
2925          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2926       ENDIF
2927!
2928!--    Organic carbon (OC) compounds
2929       IF ( index_oc > 0 )  THEN
2930          vc = 2
2931          str = ( index_oc-1 ) * nbins_aerosol + 1
2932          endi = index_oc * nbins_aerosol
2933          ic = 1
2934          DO ss = str, endi
2935             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhooc
2936             ic = ic+1
2937          ENDDO
2938          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2939       ENDIF
2940!
2941!--    Black carbon (BC)
2942       IF ( index_bc > 0 )  THEN
2943          vc = 3
2944          str = ( index_bc-1 ) * nbins_aerosol + 1 + end_subrange_1a
2945          endi = index_bc * nbins_aerosol
2946          ic = 1 + end_subrange_1a
2947          DO ss = str, endi
2948             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhobc
2949             ic = ic+1
2950          ENDDO
2951          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2952       ENDIF
2953!
2954!--    Dust (DU)
2955       IF ( index_du > 0 )  THEN
2956          vc = 4
2957          str = ( index_du-1 ) * nbins_aerosol + 1 + end_subrange_1a
2958          endi = index_du * nbins_aerosol
2959          ic = 1 + end_subrange_1a
2960          DO ss = str, endi
2961             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhodu
2962             ic = ic+1
2963          ENDDO
2964          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2965       ENDIF
2966!
2967!--    Sea salt (SS)
2968       IF ( index_ss > 0 )  THEN
2969          vc = 5
2970          str = ( index_ss-1 ) * nbins_aerosol + 1 + end_subrange_1a
2971          endi = index_ss * nbins_aerosol
2972          ic = 1 + end_subrange_1a
2973          DO ss = str, endi
2974             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoss
2975             ic = ic+1
2976          ENDDO
2977          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2978       ENDIF
2979!
2980!--    Nitrate (NO(3-)) or nitric acid HNO3
2981       IF ( index_no > 0 )  THEN
2982          vc = 6
2983          str = ( index_no-1 ) * nbins_aerosol + 1 
2984          endi = index_no * nbins_aerosol
2985          ic = 1
2986          DO ss = str, endi
2987             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhohno3
2988             ic = ic+1
2989          ENDDO
2990          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
2991       ENDIF
2992!
2993!--    Ammonium (NH(4+)) or ammonia NH3
2994       IF ( index_nh > 0 )  THEN
2995          vc = 7
2996          str = ( index_nh-1 ) * nbins_aerosol + 1
2997          endi = index_nh * nbins_aerosol
2998          ic = 1
2999          DO ss = str, endi
3000             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhonh3
3001             ic = ic+1
3002          ENDDO
3003          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
3004       ENDIF
3005!
3006!--    Water (always used)
3007       nc_h2o = get_index( prtcl,'H2O' )
3008       vc = 8
3009       str = ( nc_h2o-1 ) * nbins_aerosol + 1
3010       endi = nc_h2o * nbins_aerosol
3011       ic = 1
3012       IF ( advect_particle_water )  THEN
3013          DO ss = str, endi
3014             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2o
3015             ic = ic+1
3016          ENDDO
3017       ELSE
3018         lo_aero(1:nbins_aerosol)%volc(vc) = mclim
3019       ENDIF
3020       aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
3021!
3022!--    Number concentrations (numc) and particle sizes
3023!--    (dwet = wet diameter, core = dry volume)
3024       DO  ib = 1, nbins_aerosol
3025          lo_aero(ib)%numc = aerosol_number(ib)%conc(k,j,i)
3026          aero_old(ib)%numc = lo_aero(ib)%numc
3027          IF ( lo_aero(ib)%numc > nclim )  THEN
3028             lo_aero(ib)%dwet = ( SUM( lo_aero(ib)%volc(:) ) / lo_aero(ib)%numc / api6 )**0.33333333_wp
3029             lo_aero(ib)%core = SUM( lo_aero(ib)%volc(1:7) ) / lo_aero(ib)%numc
3030          ELSE
3031             lo_aero(ib)%dwet = lo_aero(ib)%dmid
3032             lo_aero(ib)%core = api6 * ( lo_aero(ib)%dwet )**3
3033          ENDIF
3034       ENDDO
3035!
3036!--    Calculate the ambient sizes of particles by equilibrating soluble fraction of particles with
3037!--    water using the ZSR method.
3038       in_rh = in_cw(k) / in_cs(k)
3039       IF ( prunmode==1  .OR.  .NOT. advect_particle_water )  THEN
3040          CALL equilibration( in_rh, in_t(k), lo_aero, .TRUE. )
3041       ENDIF
3042!
3043!--    Gaseous tracer concentrations in #/m3
3044       IF ( salsa_gases_from_chem )  THEN
3045!
3046!--       Convert concentrations in ppm to #/m3
3047          zgso4  = chem_species(gas_index_chem(1))%conc(k,j,i) * ppm_to_nconc(k)
3048          zghno3 = chem_species(gas_index_chem(2))%conc(k,j,i) * ppm_to_nconc(k)
3049          zgnh3  = chem_species(gas_index_chem(3))%conc(k,j,i) * ppm_to_nconc(k)
3050          zgocnv = chem_species(gas_index_chem(4))%conc(k,j,i) * ppm_to_nconc(k)
3051          zgocsv = chem_species(gas_index_chem(5))%conc(k,j,i) * ppm_to_nconc(k)
3052       ELSE
3053          zgso4  = salsa_gas(1)%conc(k,j,i)
3054          zghno3 = salsa_gas(2)%conc(k,j,i)
3055          zgnh3  = salsa_gas(3)%conc(k,j,i)
3056          zgocnv = salsa_gas(4)%conc(k,j,i)
3057          zgocsv = salsa_gas(5)%conc(k,j,i)
3058       ENDIF
3059!
3060!--    Calculate aerosol processes:
3061!--    *********************************************************************************************
3062!
3063!--    Coagulation
3064       IF ( lscoag )   THEN
3065          CALL coagulation( lo_aero, dt_salsa, in_t(k), in_p(k) )
3066       ENDIF
3067!
3068!--    Condensation
3069       IF ( lscnd )   THEN
3070          CALL condensation( lo_aero, zgso4, zgocnv, zgocsv,  zghno3, zgnh3, in_cw(k), in_cs(k),   &
3071                             in_t(k), in_p(k), dt_salsa, prtcl )
3072       ENDIF
3073!
3074!--    Deposition
3075       IF ( lsdepo )  THEN
3076          CALL deposition( lo_aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:),&
3077                           vd(k,:) )
3078       ENDIF
3079!
3080!--    Size distribution bin update
3081       IF ( lsdistupdate )   THEN
3082          CALL distr_update( lo_aero )
3083       ENDIF
3084!--    *********************************************************************************************
3085
3086       IF ( lsdepo ) sedim_vd(k,j,i,:) = vd(k,:)
3087!
3088!--    Calculate changes in concentrations
3089       DO ib = 1, nbins_aerosol
3090          aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc -   &
3091                                           aero_old(ib)%numc ) * flag
3092       ENDDO
3093
3094       IF ( index_so4 > 0 )  THEN
3095          vc = 1
3096          str = ( index_so4-1 ) * nbins_aerosol + 1
3097          endi = index_so4 * nbins_aerosol
3098          ic = 1
3099          DO ss = str, endi
3100             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3101                                            aero_old(ic)%volc(vc) ) * arhoh2so4 * flag
3102             ic = ic+1
3103          ENDDO
3104       ENDIF
3105
3106       IF ( index_oc > 0 )  THEN
3107          vc = 2
3108          str = ( index_oc-1 ) * nbins_aerosol + 1
3109          endi = index_oc * nbins_aerosol
3110          ic = 1
3111          DO ss = str, endi
3112             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3113                                            aero_old(ic)%volc(vc) ) * arhooc * flag
3114             ic = ic+1
3115          ENDDO
3116       ENDIF
3117
3118       IF ( index_bc > 0 )  THEN
3119          vc = 3
3120          str = ( index_bc-1 ) * nbins_aerosol + 1 + end_subrange_1a
3121          endi = index_bc * nbins_aerosol
3122          ic = 1 + end_subrange_1a
3123          DO ss = str, endi
3124             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3125                                            aero_old(ic)%volc(vc) ) * arhobc * flag
3126             ic = ic+1
3127          ENDDO
3128       ENDIF
3129
3130       IF ( index_du > 0 )  THEN
3131          vc = 4
3132          str = ( index_du-1 ) * nbins_aerosol + 1 + end_subrange_1a
3133          endi = index_du * nbins_aerosol
3134          ic = 1 + end_subrange_1a
3135          DO ss = str, endi
3136             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3137                                            aero_old(ic)%volc(vc) ) * arhodu * flag
3138             ic = ic+1
3139          ENDDO
3140       ENDIF
3141
3142       IF ( index_ss > 0 )  THEN
3143          vc = 5
3144          str = ( index_ss-1 ) * nbins_aerosol + 1 + end_subrange_1a
3145          endi = index_ss * nbins_aerosol
3146          ic = 1 + end_subrange_1a
3147          DO ss = str, endi
3148             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3149                                            aero_old(ic)%volc(vc) ) * arhoss * flag
3150             ic = ic+1
3151          ENDDO
3152       ENDIF
3153
3154       IF ( index_no > 0 )  THEN
3155          vc = 6
3156          str = ( index_no-1 ) * nbins_aerosol + 1
3157          endi = index_no * nbins_aerosol
3158          ic = 1
3159          DO ss = str, endi
3160             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3161                                            aero_old(ic)%volc(vc) ) * arhohno3 * flag
3162             ic = ic+1
3163          ENDDO
3164       ENDIF
3165
3166       IF ( index_nh > 0 )  THEN
3167          vc = 7
3168          str = ( index_nh-1 ) * nbins_aerosol + 1
3169          endi = index_nh * nbins_aerosol
3170          ic = 1
3171          DO ss = str, endi
3172             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3173                                            aero_old(ic)%volc(vc) ) * arhonh3 * flag
3174             ic = ic+1
3175          ENDDO
3176       ENDIF
3177
3178       IF ( advect_particle_water )  THEN
3179          nc_h2o = get_index( prtcl,'H2O' )
3180          vc = 8
3181          str = ( nc_h2o-1 ) * nbins_aerosol + 1
3182          endi = nc_h2o * nbins_aerosol
3183          ic = 1
3184          DO ss = str, endi
3185             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
3186                                            aero_old(ic)%volc(vc) ) * arhoh2o * flag
3187             ic = ic+1
3188          ENDDO
3189       ENDIF
3190       IF ( prunmode == 1 )  THEN
3191          nc_h2o = get_index( prtcl,'H2O' )
3192          vc = 8
3193          str = ( nc_h2o-1 ) * nbins_aerosol + 1
3194          endi = nc_h2o * nbins_aerosol
3195          ic = 1
3196          DO ss = str, endi
3197             aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
3198                                             aero_old(ic)%volc(vc) ) * arhoh2o )
3199             IF ( k == nzb+1 )  THEN
3200                aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
3201             ELSEIF ( k == nzt  )  THEN
3202                aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
3203                aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
3204             ENDIF
3205             ic = ic+1
3206          ENDDO
3207       ENDIF
3208!
3209!--    Condensation of precursor gases
3210       IF ( lscndgas )  THEN
3211          IF ( salsa_gases_from_chem )  THEN
3212!
3213!--          SO4 (or H2SO4)
3214             ig = gas_index_chem(1)
3215             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgso4 /               &
3216                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3217!
3218!--          HNO3
3219             ig = gas_index_chem(2)
3220             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zghno3 /              &
3221                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3222!
3223!--          NH3
3224             ig = gas_index_chem(3)
3225             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgnh3 /               &
3226                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3227!
3228!--          non-volatile OC
3229             ig = gas_index_chem(4)
3230             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgocnv /              &
3231                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3232!
3233!--          semi-volatile OC
3234             ig = gas_index_chem(5)
3235             chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgocsv /              &
3236                                            ppm_to_nconc(k) - chem_species(ig)%conc(k,j,i) ) * flag
3237
3238          ELSE
3239!
3240!--          SO4 (or H2SO4)
3241             salsa_gas(1)%conc(k,j,i) = salsa_gas(1)%conc(k,j,i) + ( zgso4 -                       &
3242                                        salsa_gas(1)%conc(k,j,i) ) * flag
3243!
3244!--          HNO3
3245             salsa_gas(2)%conc(k,j,i) = salsa_gas(2)%conc(k,j,i) + ( zghno3 -                      &
3246                                        salsa_gas(2)%conc(k,j,i) ) * flag
3247!
3248!--          NH3
3249             salsa_gas(3)%conc(k,j,i) = salsa_gas(3)%conc(k,j,i) + ( zgnh3 -                       &
3250                                        salsa_gas(3)%conc(k,j,i) ) * flag
3251!
3252!--          non-volatile OC
3253             salsa_gas(4)%conc(k,j,i) = salsa_gas(4)%conc(k,j,i) + ( zgocnv -                      &
3254                                        salsa_gas(4)%conc(k,j,i) ) * flag
3255!
3256!--          semi-volatile OC
3257             salsa_gas(5)%conc(k,j,i) = salsa_gas(5)%conc(k,j,i) + ( zgocsv -                      &
3258                                        salsa_gas(5)%conc(k,j,i) ) * flag
3259          ENDIF
3260       ENDIF
3261!
3262!--    Tendency of water vapour mixing ratio is obtained from the change in RH during SALSA run.
3263!--    This releases heat and changes pt. Assumes no temperature change during SALSA run.
3264!--    q = r / (1+r), Euler method for integration
3265!
3266       IF ( feedback_to_palm )  THEN
3267          q_p(k,j,i) = q_p(k,j,i) + 1.0_wp / ( in_cw(k) * in_adn(k) + 1.0_wp )**2 *                &
3268                       ( in_cw(k) - cw_old ) * in_adn(k) * flag
3269          pt_p(k,j,i) = pt_p(k,j,i) + alv / c_p * ( in_cw(k) - cw_old ) * in_adn(k) / ( in_cw(k) / &
3270                        in_adn(k) + 1.0_wp )**2 * pt_p(k,j,i) / in_t(k) * flag
3271       ENDIF
3272
3273    ENDDO   ! k
3274
3275!
3276!-- Set surfaces and wall fluxes due to deposition
3277    IF ( lsdepo  .AND.  lsdepo_surf  .AND.  prunmode == 3 )  THEN
3278       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
3279          CALL depo_surf( i, j, surf_def_h(0), vd, schmidt_num, kvis, in_u, .TRUE. )
3280          DO  l = 0, 3
3281             CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE. )
3282          ENDDO
3283       ELSE
3284          CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE., usm_to_depo_h )
3285          DO  l = 0, 3
3286             CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
3287                             usm_to_depo_v(l) )
3288          ENDDO
3289          CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE., lsm_to_depo_h )
3290          DO  l = 0, 3
3291             CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
3292                             lsm_to_depo_v(l) )
3293          ENDDO
3294       ENDIF
3295    ENDIF
3296
3297    IF ( prunmode < 3 )  THEN
3298       !$OMP MASTER
3299       aero = lo_aero
3300       !$OMP END MASTER
3301    END IF
3302
3303 END SUBROUTINE salsa_driver
3304
3305!------------------------------------------------------------------------------!
3306! Description:
3307! ------------
3308!> Set logical switches according to the salsa_parameters options.
3309!> Juha Tonttila, FMI, 2014
3310!> Only aerosol processes included, Mona Kurppa, UHel, 2017
3311!------------------------------------------------------------------------------!
3312 SUBROUTINE set_salsa_runtime( prunmode )
3313
3314    IMPLICIT NONE
3315
3316    INTEGER(iwp), INTENT(in) ::  prunmode
3317
3318    SELECT CASE(prunmode)
3319
3320       CASE(1) !< Initialization
3321          lscoag       = .FALSE.
3322          lscnd        = .FALSE.
3323          lscndgas     = .FALSE.
3324          lscndh2oae   = .FALSE.
3325          lsdepo       = .FALSE.
3326          lsdepo_pcm   = .FALSE.
3327          lsdepo_surf  = .FALSE.
3328          lsdistupdate = .TRUE.
3329          lspartition  = .FALSE.
3330
3331       CASE(2)  !< Spinup period
3332          lscoag      = ( .FALSE. .AND. nlcoag   )
3333          lscnd       = ( .TRUE.  .AND. nlcnd    )
3334          lscndgas    = ( .TRUE.  .AND. nlcndgas )
3335          lscndh2oae  = ( .TRUE.  .AND. nlcndh2oae )
3336
3337       CASE(3)  !< Run
3338          lscoag       = nlcoag
3339          lscnd        = nlcnd
3340          lscndgas     = nlcndgas
3341          lscndh2oae   = nlcndh2oae
3342          lsdepo       = nldepo
3343          lsdepo_pcm   = nldepo_pcm
3344          lsdepo_surf  = nldepo_surf
3345          lsdistupdate = nldistupdate
3346    END SELECT
3347
3348
3349 END SUBROUTINE set_salsa_runtime
3350 
3351!------------------------------------------------------------------------------!
3352! Description:
3353! ------------
3354!> Calculates the absolute temperature (using hydrostatic pressure), saturation
3355!> vapour pressure and mixing ratio over water, relative humidity and air
3356!> density needed in the SALSA model.
3357!> NOTE, no saturation adjustment takes place -> the resulting water vapour
3358!> mixing ratio can be supersaturated, allowing the microphysical calculations
3359!> in SALSA.
3360!
3361!> Juha Tonttila, FMI, 2014 (original SALSAthrm)
3362!> Mona Kurppa, UHel, 2017 (adjustment for PALM and only aerosol processes)
3363!------------------------------------------------------------------------------!
3364 SUBROUTINE salsa_thrm_ij( i, j, p_ij, temp_ij, cw_ij, cs_ij, adn_ij )
3365
3366    USE arrays_3d,                                                                                 &
3367        ONLY: pt, q, zu
3368
3369    USE basic_constants_and_equations_mod,                                                         &
3370        ONLY:  barometric_formula, exner_function, ideal_gas_law_rho, magnus
3371
3372    IMPLICIT NONE
3373
3374    INTEGER(iwp), INTENT(in) ::  i  !<
3375    INTEGER(iwp), INTENT(in) ::  j  !<
3376
3377    REAL(wp) ::  t_surface  !< absolute surface temperature (K)
3378
3379    REAL(wp), DIMENSION(nzb:nzt+1) ::  e_s  !< saturation vapour pressure over water (Pa)
3380
3381    REAL(wp), DIMENSION(:), INTENT(inout) ::  adn_ij   !< air density (kg/m3)
3382    REAL(wp), DIMENSION(:), INTENT(inout) ::  p_ij     !< air pressure (Pa)
3383    REAL(wp), DIMENSION(:), INTENT(inout) ::  temp_ij  !< air temperature (K)
3384
3385    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cw_ij  !< water vapour concentration (kg/m3)
3386    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cs_ij  !< saturation water vap. conc.(kg/m3)
3387!
3388!-- Pressure p_ijk (Pa) = hydrostatic pressure
3389    t_surface = pt_surface * exner_function( surface_pressure * 100.0_wp )
3390    p_ij(:) = barometric_formula( zu, t_surface, surface_pressure * 100.0_wp )
3391!
3392!-- Absolute ambient temperature (K)
3393    temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) )
3394!
3395!-- Air density
3396    adn_ij(:) = ideal_gas_law_rho( p_ij(:), temp_ij(:) )
3397!
3398!-- Water vapour concentration r_v (kg/m3)
3399    IF ( PRESENT( cw_ij ) )  THEN
3400       cw_ij(:) = ( q(:,j,i) / ( 1.0_wp - q(:,j,i) ) ) * adn_ij(:)
3401    ENDIF
3402!
3403!-- Saturation mixing ratio r_s (kg/kg) from vapour pressure at temp (Pa)
3404    IF ( PRESENT( cs_ij ) )  THEN
3405       e_s(:) = 611.0_wp * EXP( alv_d_rv * ( 3.6609E-3_wp - 1.0_wp /           &
3406                temp_ij(:) ) )! magnus( temp_ij(:) )
3407       cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:) - e_s(:) ) ) * adn_ij(:)
3408    ENDIF
3409
3410 END SUBROUTINE salsa_thrm_ij
3411
3412!------------------------------------------------------------------------------!
3413! Description:
3414! ------------
3415!> Calculates ambient sizes of particles by equilibrating soluble fraction of
3416!> particles with water using the ZSR method (Stokes and Robinson, 1966).
3417!> Method:
3418!> Following chemical components are assumed water-soluble
3419!> - (ammonium) sulphate (100%)
3420!> - sea salt (100 %)
3421!> - organic carbon (epsoc * 100%)
3422!> Exact thermodynamic considerations neglected.
3423!> - If particles contain no sea salt, calculation according to sulphate
3424!>   properties
3425!> - If contain sea salt but no sulphate, calculation according to sea salt
3426!>   properties
3427!> - If contain both sulphate and sea salt -> the molar fraction of these
3428!>   compounds determines which one of them is used as the basis of calculation.
3429!> If sulphate and sea salt coexist in a particle, it is assumed that the Cl is
3430!> replaced by sulphate; thus only either sulphate + organics or sea salt +
3431!> organics is included in the calculation of soluble fraction.
3432!> Molality parameterizations taken from Table 1 of Tang: Thermodynamic and
3433!> optical properties of mixed-salt aerosols of atmospheric importance,
3434!> J. Geophys. Res., 102 (D2), 1883-1893 (1997)
3435!
3436!> Coded by:
3437!> Hannele Korhonen (FMI) 2005
3438!> Harri Kokkola (FMI) 2006
3439!> Matti Niskanen(FMI) 2012
3440!> Anton Laakso  (FMI) 2013
3441!> Modified for the new aerosol datatype, Juha Tonttila (FMI) 2014
3442!
3443!> fxm: should sea salt form a solid particle when prh is very low (even though
3444!> it could be mixed with e.g. sulphate)?
3445!> fxm: crashes if no sulphate or sea salt
3446!> fxm: do we really need to consider Kelvin effect for subrange 2
3447!------------------------------------------------------------------------------!
3448 SUBROUTINE equilibration( prh, ptemp, paero, init )
3449
3450    IMPLICIT NONE
3451
3452    INTEGER(iwp) :: ib      !< loop index
3453    INTEGER(iwp) :: counti  !< loop index
3454
3455    LOGICAL, INTENT(in) ::  init   !< TRUE: Initialization, FALSE: Normal runtime: update water
3456                                   !< content only for 1a
3457
3458    REAL(wp) ::  zaw      !< water activity [0-1]
3459    REAL(wp) ::  zcore    !< Volume of dry particle
3460    REAL(wp) ::  zdold    !< Old diameter
3461    REAL(wp) ::  zdwet    !< Wet diameter or mean droplet diameter
3462    REAL(wp) ::  zke      !< Kelvin term in the Köhler equation
3463    REAL(wp) ::  zlwc     !< liquid water content [kg/m3-air]
3464    REAL(wp) ::  zrh      !< Relative humidity
3465
3466    REAL(wp), DIMENSION(maxspec) ::  zbinmol  !< binary molality of each components (mol/kg)
3467    REAL(wp), DIMENSION(maxspec) ::  zvpart   !< volume of chem. compounds in one particle
3468
3469    REAL(wp), INTENT(in) ::  prh    !< relative humidity [0-1]
3470    REAL(wp), INTENT(in) ::  ptemp  !< temperature (K)
3471
3472    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< aerosol properties
3473
3474    zaw       = 0.0_wp
3475    zlwc      = 0.0_wp
3476!
3477!-- Relative humidity:
3478    zrh = prh
3479    zrh = MAX( zrh, 0.05_wp )
3480    zrh = MIN( zrh, 0.98_wp)
3481!
3482!-- 1) Regime 1: sulphate and partly water-soluble OC. Done for every CALL
3483    DO  ib = start_subrange_1a, end_subrange_1a   ! size bin
3484
3485       zbinmol = 0.0_wp
3486       zdold   = 1.0_wp
3487       zke     = 1.02_wp
3488
3489       IF ( paero(ib)%numc > nclim )  THEN
3490!
3491!--       Volume in one particle
3492          zvpart = 0.0_wp
3493          zvpart(1:2) = paero(ib)%volc(1:2) / paero(ib)%numc
3494          zvpart(6:7) = paero(ib)%volc(6:7) / paero(ib)%numc
3495!
3496!--       Total volume and wet diameter of one dry particle
3497          zcore = SUM( zvpart(1:2) )
3498          zdwet = paero(ib)%dwet
3499
3500          counti = 0
3501          DO  WHILE ( ABS( zdwet / zdold - 1.0_wp ) > 1.0E-2_wp )
3502
3503             zdold = MAX( zdwet, 1.0E-20_wp )
3504             zaw = MAX( 1.0E-3_wp, zrh / zke ) ! To avoid underflow
3505!
3506!--          Binary molalities (mol/kg):
3507!--          Sulphate
3508             zbinmol(1) = 1.1065495E+2_wp - 3.6759197E+2_wp * zaw + 5.0462934E+2_wp * zaw**2 -     &
3509                          3.1543839E+2_wp * zaw**3 + 6.770824E+1_wp  * zaw**4
3510!--          Organic carbon
3511             zbinmol(2) = 1.0_wp / ( zaw * amh2o ) - 1.0_wp / amh2o
3512!--          Nitric acid
3513             zbinmol(6) = 2.306844303E+1_wp - 3.563608869E+1_wp * zaw - 6.210577919E+1_wp * zaw**2 &
3514                          + 5.510176187E+2_wp * zaw**3 - 1.460055286E+3_wp * zaw**4                &
3515                          + 1.894467542E+3_wp * zaw**5 - 1.220611402E+3_wp * zaw**6                &
3516                          + 3.098597737E+2_wp * zaw**7
3517!
3518!--          Calculate the liquid water content (kg/m3-air) using ZSR (see e.g. Eq. 10.98 in
3519!--          Seinfeld and Pandis (2006))
3520             zlwc = ( paero(ib)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / zbinmol(1) +                 &
3521                    epsoc * paero(ib)%volc(2) * ( arhooc / amoc ) / zbinmol(2) +                   &
3522                    ( paero(ib)%volc(6) * ( arhohno3/amhno3 ) ) / zbinmol(6)
3523!
3524!--          Particle wet diameter (m)
3525             zdwet = ( zlwc / paero(ib)%numc / arhoh2o / api6 + ( SUM( zvpart(6:7) ) / api6 ) +    &
3526                       zcore / api6 )**0.33333333_wp
3527!
3528!--          Kelvin effect (Eq. 10.85 in in Seinfeld and Pandis (2006)). Avoid
3529!--          overflow.
3530             zke = EXP( MIN( 50.0_wp, 4.0_wp * surfw0 * amvh2so4 / ( abo * ptemp *  zdwet ) ) )
3531
3532             counti = counti + 1
3533             IF ( counti > 1000 )  THEN
3534                message_string = 'Subrange 1: no convergence!'
3535                CALL message( 'salsa_mod: equilibration', 'PA0617', 1, 2, 0, 6, 0 )
3536             ENDIF
3537          ENDDO
3538!
3539!--       Instead of lwc, use the volume concentration of water from now on
3540!--       (easy to convert...)
3541          paero(ib)%volc(8) = zlwc / arhoh2o
3542!
3543!--       If this is initialization, update the core and wet diameter
3544          IF ( init )  THEN
3545             paero(ib)%dwet = zdwet
3546             paero(ib)%core = zcore
3547          ENDIF
3548
3549       ELSE
3550!--       If initialization
3551!--       1.2) empty bins given bin average values
3552          IF ( init )  THEN
3553             paero(ib)%dwet = paero(ib)%dmid
3554             paero(ib)%core = api6 * paero(ib)%dmid**3
3555          ENDIF
3556
3557       ENDIF
3558
3559    ENDDO  ! ib
3560!
3561!-- 2) Regime 2a: sulphate, OC, BC and sea salt
3562!--    This is done only for initialization call, otherwise the water contents
3563!--    are computed via condensation
3564    IF ( init )  THEN
3565       DO  ib = start_subrange_2a, end_subrange_2b
3566!
3567!--       Initialize
3568          zke     = 1.02_wp
3569          zbinmol = 0.0_wp
3570          zdold   = 1.0_wp
3571!
3572!--       1) Particle properties calculated for non-empty bins
3573          IF ( paero(ib)%numc > nclim )  THEN
3574!
3575!--          Volume in one particle [fxm]
3576             zvpart = 0.0_wp
3577             zvpart(1:7) = paero(ib)%volc(1:7) / paero(ib)%numc
3578!
3579!--          Total volume and wet diameter of one dry particle [fxm]
3580             zcore = SUM( zvpart(1:5) )
3581             zdwet = paero(ib)%dwet
3582
3583             counti = 0
3584             DO  WHILE ( ABS( zdwet / zdold - 1.0_wp ) > 1.0E-12_wp )
3585
3586                zdold = MAX( zdwet, 1.0E-20_wp )
3587                zaw = zrh / zke
3588!
3589!--             Binary molalities (mol/kg):
3590!--             Sulphate
3591                zbinmol(1) = 1.1065495E+2_wp - 3.6759197E+2_wp * zaw + 5.0462934E+2_wp * zaw**2 -  &
3592                             3.1543839E+2_wp * zaw**3 + 6.770824E+1_wp  * zaw**4
3593!--             Organic carbon
3594                zbinmol(2) = 1.0_wp / ( zaw * amh2o ) - 1.0_wp / amh2o
3595!--             Nitric acid
3596                zbinmol(6) = 2.306844303E+1_wp          - 3.563608869E+1_wp * zaw -                &
3597                             6.210577919E+1_wp * zaw**2 + 5.510176187E+2_wp * zaw**3 -             &
3598                             1.460055286E+3_wp * zaw**4 + 1.894467542E+3_wp * zaw**5 -             &
3599                             1.220611402E+3_wp * zaw**6 + 3.098597737E+2_wp * zaw**7 
3600!--             Sea salt (natrium chloride)
3601                zbinmol(5) = 5.875248E+1_wp - 1.8781997E+2_wp * zaw + 2.7211377E+2_wp * zaw**2 -   &
3602                             1.8458287E+2_wp * zaw**3 + 4.153689E+1_wp  * zaw**4
3603!
3604!--             Calculate the liquid water content (kg/m3-air)
3605                zlwc = ( paero(ib)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / zbinmol(1) +              &
3606                       epsoc * ( paero(ib)%volc(2) * ( arhooc / amoc ) ) / zbinmol(2) +            &
3607                       ( paero(ib)%volc(6) * ( arhohno3 / amhno3 ) ) / zbinmol(6) +                &
3608                       ( paero(ib)%volc(5) * ( arhoss / amss ) ) / zbinmol(5)
3609
3610!--             Particle wet radius (m)
3611                zdwet = ( zlwc / paero(ib)%numc / arhoh2o / api6 + ( SUM( zvpart(6:7) ) / api6 )  + &
3612                           zcore / api6 )**0.33333333_wp
3613!
3614!--             Kelvin effect (Eq. 10.85 in Seinfeld and Pandis (2006))
3615                zke = EXP( MIN( 50.0_wp, 4.0_wp * surfw0 * amvh2so4 / ( abo * zdwet * ptemp ) ) )
3616
3617                counti = counti + 1
3618                IF ( counti > 1000 )  THEN
3619                   message_string = 'Subrange 2: no convergence!'
3620                CALL message( 'salsa_mod: equilibration', 'PA0618', 1, 2, 0, 6, 0 )
3621                ENDIF
3622             ENDDO
3623!
3624!--          Liquid water content; instead of LWC use the volume concentration
3625             paero(ib)%volc(8) = zlwc / arhoh2o
3626             paero(ib)%dwet    = zdwet
3627             paero(ib)%core    = zcore
3628
3629          ELSE
3630!--          2.2) empty bins given bin average values
3631             paero(ib)%dwet = paero(ib)%dmid
3632             paero(ib)%core = api6 * paero(ib)%dmid**3
3633          ENDIF
3634
3635       ENDDO   ! ib
3636    ENDIF
3637
3638 END SUBROUTINE equilibration
3639
3640!------------------------------------------------------------------------------!
3641!> Description:
3642!> ------------
3643!> Calculation of the settling velocity vc (m/s) per aerosol size bin and
3644!> deposition on plant canopy (lsdepo_pcm).
3645!
3646!> Deposition is based on either the scheme presented in:
3647!> Zhang et al. (2001), Atmos. Environ. 35, 549-560 (includes collection due to
3648!> Brownian diffusion, impaction, interception and sedimentation; hereafter ZO1)
3649!> OR
3650!> Petroff & Zhang (2010), Geosci. Model Dev. 3, 753-769 (includes also
3651!> collection due to turbulent impaction, hereafter P10)
3652!
3653!> Equation numbers refer to equation in Jacobson (2005): Fundamentals of
3654!> Atmospheric Modeling, 2nd Edition.
3655!
3656!> Subroutine follows closely sedim_SALSA in UCLALES-SALSA written by Juha
3657!> Tonttila (KIT/FMI) and Zubair Maalick (UEF).
3658!> Rewritten to PALM by Mona Kurppa (UH), 2017.
3659!
3660!> Call for grid point i,j,k
3661!------------------------------------------------------------------------------!
3662
3663 SUBROUTINE deposition( paero, tk, adn, mag_u, lad, kvis, schmidt_num, vc )
3664
3665    USE plant_canopy_model_mod,                                                                    &
3666        ONLY:  cdc
3667
3668    IMPLICIT NONE
3669
3670    INTEGER(iwp) ::  ib   !< loop index
3671    INTEGER(iwp) ::  ic   !< loop index
3672
3673    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
3674    REAL(wp) ::  avis              !< molecular viscocity of air (kg/(m*s))
3675    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
3676    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
3677    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
3678    REAL(wp) ::  c_interception    !< coefficient for interception
3679    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
3680    REAL(wp) ::  depo              !< deposition velocity (m/s)
3681    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
3682    REAL(wp) ::  lambda            !< molecular mean free path (m)
3683    REAL(wp) ::  mdiff             !< particle diffusivity coefficient
3684    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
3685                                   !< Table 3 in Z01
3686    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
3687    REAL(wp) ::  pdn               !< particle density (kg/m3)
3688    REAL(wp) ::  ustar             !< friction velocity (m/s)
3689    REAL(wp) ::  va                !< thermal speed of an air molecule (m/s)
3690
3691    REAL(wp), INTENT(in) ::  adn    !< air density (kg/m3)
3692    REAL(wp), INTENT(in) ::  lad    !< leaf area density (m2/m3)
3693    REAL(wp), INTENT(in) ::  mag_u  !< wind velocity (m/s)
3694    REAL(wp), INTENT(in) ::  tk     !< abs.temperature (K)
3695
3696    REAL(wp), INTENT(inout) ::  kvis   !< kinematic viscosity of air (m2/s)
3697
3698    REAL(wp), DIMENSION(nbins_aerosol) ::  beta   !< Cunningham slip-flow correction factor
3699    REAL(wp), DIMENSION(nbins_aerosol) ::  Kn     !< Knudsen number
3700    REAL(wp), DIMENSION(nbins_aerosol) ::  zdwet  !< wet diameter (m)
3701
3702    REAL(wp), DIMENSION(:), INTENT(inout) ::  schmidt_num  !< particle Schmidt number
3703    REAL(wp), DIMENSION(:), INTENT(inout) ::  vc  !< critical fall speed i.e. settling velocity of
3704                                                  !< an aerosol particle (m/s)
3705
3706    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< aerosol properties
3707!
3708!-- Initialise
3709    depo  = 0.0_wp
3710    pdn   = 1500.0_wp    ! default value
3711    ustar = 0.0_wp
3712!
3713!-- Molecular viscosity of air (Eq. 4.54)
3714    avis = 1.8325E-5_wp * ( 416.16_wp / ( tk + 120.0_wp ) ) * ( tk / 296.16_wp )**1.5_wp
3715!
3716!-- Kinematic viscosity (Eq. 4.55)
3717    kvis =  avis / adn
3718!
3719!-- Thermal velocity of an air molecule (Eq. 15.32)
3720    va = SQRT( 8.0_wp * abo * tk / ( pi * am_airmol ) )
3721!
3722!-- Mean free path (m) (Eq. 15.24)
3723    lambda = 2.0_wp * avis / ( adn * va )
3724!
3725!-- Particle wet diameter (m)
3726    zdwet = paero(:)%dwet
3727!
3728!-- Knudsen number (Eq. 15.23)
3729    Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
3730!
3731!-- Cunningham slip-flow correction (Eq. 15.30)
3732    beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
3733!
3734!-- Critical fall speed i.e. settling velocity  (Eq. 20.4)
3735    vc = MIN( 1.0_wp, zdwet**2 * ( pdn - adn ) * g * beta / ( 18.0_wp * avis ) )
3736!
3737!-- Deposition on vegetation
3738    IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
3739!
3740!--    Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
3741       alpha   = alpha_z01(depo_pcm_type_num)
3742       gamma   = gamma_z01(depo_pcm_type_num)
3743       par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
3744!
3745!--    Deposition efficiencies from Table 1. Constants from Table 2.
3746       par_l            = l_p10(depo_pcm_type_num) * 0.01_wp
3747       c_brownian_diff  = c_b_p10(depo_pcm_type_num)
3748       c_interception   = c_in_p10(depo_pcm_type_num)
3749       c_impaction      = c_im_p10(depo_pcm_type_num)
3750       beta_im          = beta_im_p10(depo_pcm_type_num)
3751       c_turb_impaction = c_it_p10(depo_pcm_type_num)
3752
3753       DO  ib = 1, nbins_aerosol
3754
3755          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
3756
3757!--       Particle diffusivity coefficient (Eq. 15.29)
3758          mdiff = ( abo * tk * beta(ib) ) / ( 3.0_wp * pi * avis * zdwet(ib) )
3759!
3760!--       Particle Schmidt number (Eq. 15.36)
3761          schmidt_num(ib) = kvis / mdiff
3762!
3763!--       Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
3764          ustar = SQRT( cdc ) * mag_u
3765          SELECT CASE ( depo_pcm_par_num )
3766
3767             CASE ( 1 )   ! Zhang et al. (2001)
3768                CALL depo_vel_Z01( vc(ib), ustar, schmidt_num(ib), paero(ib)%dwet, alpha,  gamma,  &
3769                                   par_a, depo )
3770             CASE ( 2 )   ! Petroff & Zhang (2010)
3771                CALL depo_vel_P10( vc(ib), mag_u, ustar, kvis, schmidt_num(ib), paero(ib)%dwet,    &
3772                                   par_l, c_brownian_diff, c_interception, c_impaction, beta_im,   &
3773                                   c_turb_impaction, depo )
3774          END SELECT
3775!
3776!--       Calculate the change in concentrations
3777          paero(ib)%numc = paero(ib)%numc - depo * lad * paero(ib)%numc * dt_salsa
3778          DO  ic = 1, maxspec+1
3779             paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa
3780          ENDDO
3781       ENDDO
3782
3783    ENDIF
3784
3785 END SUBROUTINE deposition
3786
3787!------------------------------------------------------------------------------!
3788! Description:
3789! ------------
3790!> Calculate deposition velocity (m/s) based on Zhan et al. (2001, case 1).
3791!------------------------------------------------------------------------------!
3792
3793 SUBROUTINE depo_vel_Z01( vc, ustar, schmidt_num, diameter, alpha, gamma, par_a, depo )
3794
3795    IMPLICIT NONE
3796
3797    REAL(wp) ::  rs                !< overall quasi-laminar resistance for particles
3798    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
3799
3800    REAL(wp), INTENT(in) ::  alpha        !< parameter, Table 3 in Z01
3801    REAL(wp), INTENT(in) ::  gamma        !< parameter, Table 3 in Z01
3802    REAL(wp), INTENT(in) ::  par_a        !< parameter A for the characteristic diameter of
3803                                          !< collectors, Table 3 in Z01
3804    REAL(wp), INTENT(in) ::  diameter     !< particle diameter
3805    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
3806    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
3807    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
3808
3809    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
3810
3811    IF ( par_a > 0.0_wp )  THEN
3812!
3813!--    Initialise
3814       rs = 0.0_wp
3815!
3816!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
3817       stokes_num = vc * ustar / ( g * par_a )
3818!
3819!--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)
3820       rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -stokes_num**0.5_wp ) *                &
3821                 ( schmidt_num**( -gamma ) + ( stokes_num / ( alpha + stokes_num ) )**2 +          &
3822                 0.5_wp * ( diameter / par_a )**2 ) ) )
3823
3824       depo = rs + vc
3825
3826    ELSE
3827       depo = 0.0_wp
3828    ENDIF
3829
3830 END SUBROUTINE depo_vel_Z01
3831
3832!------------------------------------------------------------------------------!
3833! Description:
3834! ------------
3835!> Calculate deposition velocity (m/s) based on Petroff & Zhang (2010, case 2).
3836!------------------------------------------------------------------------------!
3837
3838 SUBROUTINE depo_vel_P10( vc, mag_u, ustar, kvis_a, schmidt_num, diameter, par_l, c_brownian_diff, &
3839                          c_interception, c_impaction, beta_im, c_turb_impaction, depo )
3840
3841    IMPLICIT NONE
3842
3843    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
3844    REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
3845    REAL(wp) ::  v_bd              !< deposition velocity due to Brownian diffusion
3846    REAL(wp) ::  v_im              !< deposition velocity due to impaction
3847    REAL(wp) ::  v_in              !< deposition velocity due to interception
3848    REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
3849
3850    REAL(wp), INTENT(in) ::  beta_im           !< parameter for turbulent impaction
3851    REAL(wp), INTENT(in) ::  c_brownian_diff   !< coefficient for Brownian diffusion
3852    REAL(wp), INTENT(in) ::  c_impaction       !< coefficient for inertial impaction
3853    REAL(wp), INTENT(in) ::  c_interception    !< coefficient for interception
3854    REAL(wp), INTENT(in) ::  c_turb_impaction  !< coefficient for turbulent impaction
3855    REAL(wp), INTENT(in) ::  kvis_a       !< kinematic viscosity of air (m2/s)
3856    REAL(wp), INTENT(in) ::  mag_u        !< wind velocity (m/s)
3857    REAL(wp), INTENT(in) ::  par_l        !< obstacle characteristic dimension in P10
3858    REAL(wp), INTENT(in) ::  diameter       !< particle diameter
3859    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
3860    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
3861    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
3862
3863    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
3864
3865    IF ( par_l > 0.0_wp )  THEN
3866!
3867!--    Initialise
3868       tau_plus = 0.0_wp
3869       v_bd     = 0.0_wp
3870       v_im     = 0.0_wp
3871       v_in     = 0.0_wp
3872       v_it     = 0.0_wp
3873!
3874!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
3875       stokes_num = vc * ustar / ( g * par_l )
3876!
3877!--    Non-dimensional relexation time of the particle on top of canopy
3878       tau_plus = vc * ustar**2 / ( kvis_a * g )
3879!
3880!--    Brownian diffusion
3881       v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                          &
3882              ( mag_u * par_l / kvis_a )**( -0.5_wp )
3883!
3884!--    Interception
3885       v_in = mag_u * c_interception * diameter / par_l *                                          &
3886              ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) )
3887!
3888!--    Impaction: Petroff (2009) Eq. 18
3889       v_im = mag_u * c_impaction * ( stokes_num / ( stokes_num + beta_im ) )**2
3890!
3891!--    Turbulent impaction
3892       IF ( tau_plus < 20.0_wp )  THEN
3893          v_it = 2.5E-3_wp * c_turb_impaction * tau_plus**2
3894       ELSE
3895          v_it = c_turb_impaction
3896       ENDIF
3897
3898       depo = ( v_bd + v_in + v_im + v_it + vc )
3899
3900    ELSE
3901       depo = 0.0_wp
3902    ENDIF
3903
3904 END SUBROUTINE depo_vel_P10
3905
3906!------------------------------------------------------------------------------!
3907! Description:
3908! ------------
3909!> Calculate the dry deposition on horizontal and vertical surfaces. Implement
3910!> as a surface flux.
3911!> @todo aerodynamic resistance ignored for now (not important for
3912!        high-resolution simulations)
3913!------------------------------------------------------------------------------!
3914 SUBROUTINE depo_surf( i, j, surf, vc, schmidt_num, kvis, mag_u, norm, match_array )
3915
3916    USE arrays_3d,                                                                                 &
3917        ONLY: rho_air_zw
3918
3919    USE surface_mod,                                                                               &
3920        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
3921
3922    IMPLICIT NONE
3923
3924    INTEGER(iwp) ::  ib      !< loop index
3925    INTEGER(iwp) ::  ic      !< loop index
3926    INTEGER(iwp) ::  icc     !< additional loop index
3927    INTEGER(iwp) ::  k       !< loop index
3928    INTEGER(iwp) ::  m       !< loop index
3929    INTEGER(iwp) ::  surf_e  !< End index of surface elements at (j,i)-gridpoint
3930    INTEGER(iwp) ::  surf_s  !< Start index of surface elements at (j,i)-gridpoint
3931
3932    INTEGER(iwp), INTENT(in) ::  i  !< loop index
3933    INTEGER(iwp), INTENT(in) ::  j  !< loop index
3934
3935    LOGICAL, INTENT(in) ::  norm   !< to normalise or not
3936
3937    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
3938    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
3939    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
3940    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
3941    REAL(wp) ::  c_interception    !< coefficient for interception
3942    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
3943    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
3944    REAL(wp) ::  norm_fac          !< normalisation factor (usually air density)
3945    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
3946                                   !< Table 3 in Z01
3947    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
3948    REAL(wp) ::  rs                !< the overall quasi-laminar resistance for particles
3949    REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
3950    REAL(wp) ::  v_bd              !< deposition velocity due to Brownian diffusion
3951    REAL(wp) ::  v_im              !< deposition velocity due to impaction
3952    REAL(wp) ::  v_in              !< deposition velocity due to interception
3953    REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
3954
3955    REAL(wp), DIMENSION(nbins_aerosol) ::  depo      !< deposition efficiency
3956    REAL(wp), DIMENSION(nbins_aerosol) ::  depo_sum  !< sum of deposition efficiencies
3957
3958    REAL(wp), DIMENSION(:), INTENT(in) ::  kvis   !< kinematic viscosity of air (m2/s)
3959    REAL(wp), DIMENSION(:), INTENT(in) ::  mag_u  !< wind velocity (m/s)
3960
3961    REAL(wp), DIMENSION(:,:), INTENT(in) ::  schmidt_num   !< particle Schmidt number
3962    REAL(wp), DIMENSION(:,:), INTENT(in) ::  vc            !< terminal velocity (m/s)
3963
3964    TYPE(match_surface), INTENT(in), OPTIONAL ::  match_array  !< match the deposition module and
3965                                                               !< LSM/USM surfaces
3966    TYPE(surf_type), INTENT(inout) :: surf                     !< respective surface type
3967!
3968!-- Initialise
3969    depo     = 0.0_wp
3970    depo_sum = 0.0_wp
3971    rs       = 0.0_wp
3972    surf_s   = surf%start_index(j,i)
3973    surf_e   = surf%end_index(j,i)
3974    tau_plus = 0.0_wp
3975    v_bd     = 0.0_wp
3976    v_im     = 0.0_wp
3977    v_in     = 0.0_wp
3978    v_it     = 0.0_wp
3979!
3980!-- Model parameters for the land use category. If LSM or USM is applied, import
3981!-- characteristics. Otherwise, apply surface type "urban".
3982    alpha   = alpha_z01(luc_urban)
3983    gamma   = gamma_z01(luc_urban)
3984    par_a   = A_z01(luc_urban, season) * 1.0E-3_wp
3985
3986    par_l            = l_p10(luc_urban) * 0.01_wp
3987    c_brownian_diff  = c_b_p10(luc_urban)
3988    c_interception   = c_in_p10(luc_urban)
3989    c_impaction      = c_im_p10(luc_urban)
3990    beta_im          = beta_im_p10(luc_urban)
3991    c_turb_impaction = c_it_p10(luc_urban)
3992
3993
3994    IF ( PRESENT( match_array ) )  THEN  ! land or urban surface model
3995
3996       DO  m = surf_s, surf_e
3997
3998          k = surf%k(m)
3999          norm_fac = 1.0_wp
4000
4001          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
4002
4003          IF ( match_array%match_lupg(m) > 0 )  THEN
4004             alpha = alpha_z01( match_array%match_lupg(m) )
4005             gamma = gamma_z01( match_array%match_lupg(m) )
4006             par_a = A_z01( match_array%match_lupg(m), season ) * 1.0E-3_wp
4007
4008             beta_im          = beta_im_p10( match_array%match_lupg(m) )
4009             c_brownian_diff  = c_b_p10( match_array%match_lupg(m) )
4010             c_impaction      = c_im_p10( match_array%match_lupg(m) )
4011             c_interception   = c_in_p10( match_array%match_lupg(m) )
4012             c_turb_impaction = c_it_p10( match_array%match_lupg(m) )
4013             par_l            = l_p10( match_array%match_lupg(m) ) * 0.01_wp
4014
4015             DO  ib = 1, nbins_aerosol
4016                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
4017                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
4018
4019                SELECT CASE ( depo_surf_par_num )
4020
4021                   CASE ( 1 )
4022                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
4023                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
4024                   CASE ( 2 )
4025                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
4026                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
4027                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
4028                                         c_turb_impaction, depo(ib) )
4029                END SELECT
4030             ENDDO
4031             depo_sum = depo_sum + surf%frac(ind_pav_green,m) * depo
4032          ENDIF
4033
4034          IF ( match_array%match_luvw(m) > 0 )  THEN
4035             alpha = alpha_z01( match_array%match_luvw(m) )
4036             gamma = gamma_z01( match_array%match_luvw(m) )
4037             par_a = A_z01( match_array%match_luvw(m), season ) * 1.0E-3_wp
4038
4039             beta_im          = beta_im_p10( match_array%match_luvw(m) )
4040             c_brownian_diff  = c_b_p10( match_array%match_luvw(m) )
4041             c_impaction      = c_im_p10( match_array%match_luvw(m) )
4042             c_interception   = c_in_p10( match_array%match_luvw(m) )
4043             c_turb_impaction = c_it_p10( match_array%match_luvw(m) )
4044             par_l            = l_p10( match_array%match_luvw(m) ) * 0.01_wp
4045
4046             DO  ib = 1, nbins_aerosol
4047                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
4048                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
4049
4050                SELECT CASE ( depo_surf_par_num )
4051
4052                   CASE ( 1 )
4053                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
4054                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
4055                   CASE ( 2 )
4056                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
4057                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
4058                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
4059                                         c_turb_impaction, depo(ib) )
4060                END SELECT
4061             ENDDO
4062             depo_sum = depo_sum + surf%frac(ind_veg_wall,m) * depo
4063          ENDIF
4064
4065          IF ( match_array%match_luww(m) > 0 )  THEN
4066             alpha = alpha_z01( match_array%match_luww(m) )
4067             gamma = gamma_z01( match_array%match_luww(m) )
4068             par_a = A_z01( match_array%match_luww(m), season ) * 1.0E-3_wp
4069
4070             beta_im          = beta_im_p10( match_array%match_luww(m) )
4071             c_brownian_diff  = c_b_p10( match_array%match_luww(m) )
4072             c_impaction      = c_im_p10( match_array%match_luww(m) )
4073             c_interception   = c_in_p10( match_array%match_luww(m) )
4074             c_turb_impaction = c_it_p10( match_array%match_luww(m) )
4075             par_l            = l_p10( match_array%match_luww(m) ) * 0.01_wp
4076
4077             DO  ib = 1, nbins_aerosol
4078                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
4079                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
4080
4081                SELECT CASE ( depo_surf_par_num )
4082
4083                   CASE ( 1 )
4084                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
4085                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
4086                   CASE ( 2 )
4087                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
4088                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
4089                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
4090                                         c_turb_impaction, depo(ib) )
4091                END SELECT
4092             ENDDO
4093             depo_sum = depo_sum + surf%frac(ind_wat_win,m) * depo
4094          ENDIF
4095
4096          DO  ib = 1, nbins_aerosol
4097             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) )  CYCLE
4098!
4099!--          Calculate changes in surface fluxes due to dry deposition
4100             IF ( include_emission )  THEN
4101                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
4102                                   depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
4103                DO  ic = 1, ncomponents_mass
4104                   icc = ( ic - 1 ) * nbins_aerosol + ib
4105                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
4106                                       depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
4107                ENDDO  ! ic
4108             ELSE
4109                surf%answs(m,ib) = -depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
4110                DO  ic = 1, ncomponents_mass
4111                   icc = ( ic - 1 ) * nbins_aerosol + ib
4112                   surf%amsws(m,icc) = -depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
4113                ENDDO  ! ic
4114             ENDIF
4115          ENDDO  ! ib
4116
4117       ENDDO
4118
4119    ELSE  ! default surfaces
4120
4121       DO  m = surf_s, surf_e
4122
4123          k = surf%k(m)
4124          norm_fac = 1.0_wp
4125
4126          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
4127
4128          DO  ib = 1, nbins_aerosol
4129             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                        &
4130                  schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
4131
4132             SELECT CASE ( depo_surf_par_num )
4133
4134                CASE ( 1 )
4135                   CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),                 &
4136                                      ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
4137                CASE ( 2 )
4138                   CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),               &
4139                                      schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,                &
4140                                      c_brownian_diff, c_interception, c_impaction, beta_im,       &
4141                                      c_turb_impaction, depo(ib) )
4142             END SELECT
4143!
4144!--          Calculate changes in surface fluxes due to dry deposition
4145             IF ( include_emission )  THEN
4146                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
4147                                   depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
4148                DO  ic = 1, ncomponents_mass
4149                   icc = ( ic - 1 ) * nbins_aerosol + ib
4150                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
4151                                       depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
4152                ENDDO  ! ic
4153             ELSE
4154                surf%answs(m,ib) = -depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
4155                DO  ic = 1, ncomponents_mass
4156                   icc = ( ic - 1 ) * nbins_aerosol + ib
4157                   surf%amsws(m,icc) = -depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
4158                ENDDO  ! ic
4159             ENDIF
4160          ENDDO  ! ib
4161       ENDDO
4162
4163    ENDIF
4164
4165 END SUBROUTINE depo_surf
4166
4167!------------------------------------------------------------------------------!
4168! Description:
4169! ------------
4170!> Calculates particle loss and change in size distribution due to (Brownian)
4171!> coagulation. Only for particles with dwet < 30 micrometres.
4172!
4173!> Method:
4174!> Semi-implicit, non-iterative method: (Jacobson, 1994)
4175!> Volume concentrations of the smaller colliding particles added to the bin of
4176!> the larger colliding particles. Start from first bin and use the updated
4177!> number and volume for calculation of following bins. NB! Our bin numbering
4178!> does not follow particle size in subrange 2.
4179!
4180!> Schematic for bin numbers in different subranges:
4181!>             1                            2
4182!>    +-------------------------------------------+
4183!>  a | 1 | 2 | 3 || 4 | 5 | 6 | 7 |  8 |  9 | 10||
4184!>  b |           ||11 |12 |13 |14 | 15 | 16 | 17||
4185!>    +-------------------------------------------+
4186!
4187!> Exact coagulation coefficients for each pressure level are scaled according
4188!> to current particle wet size (linear scaling).
4189!> Bins are organized in terms of the dry size of the condensation nucleus,
4190!> while coagulation kernell is calculated with the actual hydrometeor
4191!> size.
4192!
4193!> Called from salsa_driver
4194!> fxm: Process selection should be made smarter - now just lots of IFs inside
4195!>      loops
4196!
4197!> Coded by:
4198!> Hannele Korhonen (FMI) 2005
4199!> Harri Kokkola (FMI) 2006
4200!> Tommi Bergman (FMI) 2012
4201!> Matti Niskanen(FMI) 2012
4202!> Anton Laakso  (FMI) 2013
4203!> Juha Tonttila (FMI) 2014
4204!------------------------------------------------------------------------------!
4205 SUBROUTINE coagulation( paero, ptstep, ptemp, ppres )
4206
4207    IMPLICIT NONE
4208
4209    INTEGER(iwp) ::  index_2a !< corresponding bin in subrange 2a
4210    INTEGER(iwp) ::  index_2b !< corresponding bin in subrange 2b
4211    INTEGER(iwp) ::  ib       !< loop index
4212    INTEGER(iwp) ::  ll       !< loop index
4213    INTEGER(iwp) ::  mm       !< loop index
4214    INTEGER(iwp) ::  nn       !< loop index
4215
4216    REAL(wp) ::  pressi          !< pressure
4217    REAL(wp) ::  temppi          !< temperature
4218    REAL(wp) ::  zdpart_mm       !< diameter of particle (m)
4219    REAL(wp) ::  zdpart_nn       !< diameter of particle (m)
4220    REAL(wp) ::  zminusterm      !< coagulation loss in a bin (1/s)
4221
4222    REAL(wp), INTENT(in) ::  ppres  !< ambient pressure (Pa)
4223    REAL(wp), INTENT(in) ::  ptemp  !< ambient temperature (K)
4224    REAL(wp), INTENT(in) ::  ptstep !< time step (s)
4225
4226    REAL(wp), DIMENSION(nbins_aerosol) ::  zmpart     !< approximate mass of particles (kg)
4227    REAL(wp), DIMENSION(maxspec+1)     ::  zplusterm  !< coagulation gain in a bin (for each
4228                                                      !< chemical compound)
4229    REAL(wp), DIMENSION(nbins_aerosol,nbins_aerosol) ::  zcc  !< updated coagulation coeff. (m3/s)
4230
4231    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< Aerosol properties
4232
4233    zdpart_mm = 0.0_wp
4234    zdpart_nn = 0.0_wp
4235!
4236!-- 1) Coagulation to coarse mode calculated in a simplified way:
4237!--    CoagSink ~ Dp in continuum subrange --> 'effective' number conc. of coarse particles
4238
4239!-- 2) Updating coagulation coefficients
4240!
4241!-- Aerosol mass (kg). Density of 1500 kg/m3 assumed
4242    zmpart(1:end_subrange_2b) = api6 * ( MIN( paero(1:end_subrange_2b)%dwet, 30.0E-6_wp )**3 )     &
4243                                * 1500.0_wp
4244    temppi = ptemp
4245    pressi = ppres
4246    zcc    = 0.0_wp
4247!
4248!-- Aero-aero coagulation
4249    DO  mm = 1, end_subrange_2b   ! smaller colliding particle
4250       IF ( paero(mm)%numc < ( 2.0_wp * nclim ) )  CYCLE
4251       DO  nn = mm, end_subrange_2b   ! larger colliding particle
4252          IF ( paero(nn)%numc < ( 2.0_wp * nclim ) )  CYCLE
4253
4254          zdpart_mm = MIN( paero(mm)%dwet, 30.0E-6_wp )     ! Limit to 30 um
4255          zdpart_nn = MIN( paero(nn)%dwet, 30.0E-6_wp )     ! Limit to 30 um
4256!
4257!--       Coagulation coefficient of particles (m3/s)
4258          zcc(mm,nn) = coagc( zdpart_mm, zdpart_nn, zmpart(mm), zmpart(nn), temppi, pressi )
4259          zcc(nn,mm) = zcc(mm,nn)
4260       ENDDO
4261    ENDDO
4262
4263!
4264!-- 3) New particle and volume concentrations after coagulation:
4265!--    Calculated according to Jacobson (2005) eq. 15.9
4266!
4267!-- Aerosols in subrange 1a:
4268    DO  ib = start_subrange_1a, end_subrange_1a
4269       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
4270       zminusterm   = 0.0_wp
4271       zplusterm(:) = 0.0_wp
4272!
4273!--    Particles lost by coagulation with larger aerosols
4274       DO  ll = ib+1, end_subrange_2b
4275          zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4276       ENDDO
4277!
4278!--    Coagulation gain in a bin: change in volume conc. (cm3/cm3):
4279       DO ll = start_subrange_1a, ib - 1
4280          zplusterm(1:2) = zplusterm(1:2) + zcc(ll,ib) * paero(ll)%volc(1:2)
4281          zplusterm(6:7) = zplusterm(6:7) + zcc(ll,ib) * paero(ll)%volc(6:7)
4282          zplusterm(8)   = zplusterm(8)   + zcc(ll,ib) * paero(ll)%volc(8)
4283       ENDDO
4284!
4285!--    Volume and number concentrations after coagulation update [fxm]
4286       paero(ib)%volc(1:2) = ( paero(ib)%volc(1:2) + ptstep * zplusterm(1:2) * paero(ib)%numc ) /  &
4287                            ( 1.0_wp + ptstep * zminusterm )
4288       paero(ib)%volc(6:8) = ( paero(ib)%volc(6:8) + ptstep * zplusterm(6:8) * paero(ib)%numc ) /  &
4289                            ( 1.0_wp + ptstep * zminusterm )
4290       paero(ib)%numc = paero(ib)%numc / ( 1.0_wp + ptstep * zminusterm + 0.5_wp * ptstep *        &
4291                        zcc(ib,ib) * paero(ib)%numc )
4292    ENDDO
4293!
4294!-- Aerosols in subrange 2a:
4295    DO  ib = start_subrange_2a, end_subrange_2a
4296       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
4297       zminusterm   = 0.0_wp
4298       zplusterm(:) = 0.0_wp
4299!
4300!--    Find corresponding size bin in subrange 2b
4301       index_2b = ib - start_subrange_2a + start_subrange_2b
4302!
4303!--    Particles lost by larger particles in 2a
4304       DO  ll = ib+1, end_subrange_2a
4305          zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4306       ENDDO
4307!
4308!--    Particles lost by larger particles in 2b
4309       IF ( .NOT. no_insoluble )  THEN
4310          DO  ll = index_2b+1, end_subrange_2b
4311             zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4312          ENDDO
4313       ENDIF
4314!
4315!--    Particle volume gained from smaller particles in subranges 1, 2a and 2b
4316       DO  ll = start_subrange_1a, ib-1
4317          zplusterm(1:2) = zplusterm(1:2) + zcc(ll,ib) * paero(ll)%volc(1:2)
4318          zplusterm(6:8) = zplusterm(6:8) + zcc(ll,ib) * paero(ll)%volc(6:8)
4319       ENDDO
4320!
4321!--    Particle volume gained from smaller particles in 2a
4322!--    (Note, for components not included in the previous loop!)
4323       DO  ll = start_subrange_2a, ib-1
4324          zplusterm(3:5) = zplusterm(3:5) + zcc(ll,ib)*paero(ll)%volc(3:5)
4325       ENDDO
4326!
4327!--    Particle volume gained from smaller (and equal) particles in 2b
4328       IF ( .NOT. no_insoluble )  THEN
4329          DO  ll = start_subrange_2b, index_2b
4330             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,ib) * paero(ll)%volc(1:8)
4331          ENDDO
4332       ENDIF
4333!
4334!--    Volume and number concentrations after coagulation update [fxm]
4335       paero(ib)%volc(1:8) = ( paero(ib)%volc(1:8) + ptstep * zplusterm(1:8) * paero(ib)%numc ) /  &
4336                            ( 1.0_wp + ptstep * zminusterm )
4337       paero(ib)%numc = paero(ib)%numc / ( 1.0_wp + ptstep * zminusterm + 0.5_wp * ptstep *        &
4338                        zcc(ib,ib) * paero(ib)%numc )
4339    ENDDO
4340!
4341!-- Aerosols in subrange 2b:
4342    IF ( .NOT. no_insoluble )  THEN
4343       DO  ib = start_subrange_2b, end_subrange_2b
4344          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
4345          zminusterm   = 0.0_wp
4346          zplusterm(:) = 0.0_wp
4347!
4348!--       Find corresponding size bin in subsubrange 2a
4349          index_2a = ib - start_subrange_2b + start_subrange_2a
4350!
4351!--       Particles lost to larger particles in subranges 2b
4352          DO  ll = ib + 1, end_subrange_2b
4353             zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4354          ENDDO
4355!
4356!--       Particles lost to larger and equal particles in 2a
4357          DO  ll = index_2a, end_subrange_2a
4358             zminusterm = zminusterm + zcc(ib,ll) * paero(ll)%numc
4359          ENDDO
4360!
4361!--       Particle volume gained from smaller particles in subranges 1 & 2a
4362          DO  ll = start_subrange_1a, index_2a - 1
4363             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,ib) * paero(ll)%volc(1:8)
4364          ENDDO
4365!
4366!--       Particle volume gained from smaller particles in 2b
4367          DO  ll = start_subrange_2b, ib - 1
4368             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,ib) * paero(ll)%volc(1:8)
4369          ENDDO
4370!
4371!--       Volume and number concentrations after coagulation update [fxm]
4372          paero(ib)%volc(1:8) = ( paero(ib)%volc(1:8) + ptstep * zplusterm(1:8) * paero(ib)%numc ) &
4373                                / ( 1.0_wp + ptstep * zminusterm )
4374          paero(ib)%numc = paero(ib)%numc / ( 1.0_wp + ptstep * zminusterm + 0.5_wp * ptstep *     &
4375                           zcc(ib,ib) * paero(ib)%numc )
4376       ENDDO
4377    ENDIF
4378
4379 END SUBROUTINE coagulation
4380
4381!------------------------------------------------------------------------------!
4382! Description:
4383! ------------
4384!> Calculation of coagulation coefficients. Extended version of the function
4385!> originally found in mo_salsa_init.
4386!
4387!> J. Tonttila, FMI, 05/2014
4388!------------------------------------------------------------------------------!
4389 REAL(wp) FUNCTION coagc( diam1, diam2, mass1, mass2, temp, pres )
4390
4391    IMPLICIT NONE
4392
4393    REAL(wp) ::  fmdist  !< distance of flux matching (m)
4394    REAL(wp) ::  knud_p  !< particle Knudsen number
4395    REAL(wp) ::  mdiam   !< mean diameter of colliding particles (m)
4396    REAL(wp) ::  mfp     !< mean free path of air molecules (m)
4397    REAL(wp) ::  visc    !< viscosity of air (kg/(m s))
4398
4399    REAL(wp), INTENT(in) ::  diam1  !< diameter of colliding particle 1 (m)
4400    REAL(wp), INTENT(in) ::  diam2  !< diameter of colliding particle 2 (m)
4401    REAL(wp), INTENT(in) ::  mass1  !< mass of colliding particle 1 (kg)
4402    REAL(wp), INTENT(in) ::  mass2  !< mass of colliding particle 2 (kg)
4403    REAL(wp), INTENT(in) ::  pres   !< ambient pressure (Pa?) [fxm]
4404    REAL(wp), INTENT(in) ::  temp   !< ambient temperature (K)
4405
4406    REAL(wp), DIMENSION (2) ::  beta    !< Cunningham correction factor
4407    REAL(wp), DIMENSION (2) ::  dfpart  !< particle diffusion coefficient (m2/s)
4408    REAL(wp), DIMENSION (2) ::  diam    !< diameters of particles (m)
4409    REAL(wp), DIMENSION (2) ::  flux    !< flux in continuum and free molec. regime (m/s)
4410    REAL(wp), DIMENSION (2) ::  knud    !< particle Knudsen number
4411    REAL(wp), DIMENSION (2) ::  mpart   !< masses of particles (kg)
4412    REAL(wp), DIMENSION (2) ::  mtvel   !< particle mean thermal velocity (m/s)
4413    REAL(wp), DIMENSION (2) ::  omega   !< particle mean free path
4414    REAL(wp), DIMENSION (2) ::  tva     !< temporary variable (m)
4415!
4416!-- Initialisation
4417    coagc   = 0.0_wp
4418!
4419!-- 1) Initializing particle and ambient air variables
4420    diam  = (/ diam1, diam2 /) !< particle diameters (m)
4421    mpart = (/ mass1, mass2 /) !< particle masses (kg)
4422!
4423!-- Viscosity of air (kg/(m s))
4424    visc = ( 7.44523E-3_wp * temp ** 1.5_wp ) / ( 5093.0_wp * ( temp + 110.4_wp ) )
4425!
4426!-- Mean free path of air (m)
4427    mfp = ( 1.656E-10_wp * temp + 1.828E-8_wp ) * ( p_0 + 1325.0_wp ) / pres
4428!
4429!-- 2) Slip correction factor for small particles
4430    knud = 2.0_wp * EXP( LOG(mfp) - LOG(diam) )! Knudsen number for air (15.23)
4431!
4432!-- Cunningham correction factor (Allen and Raabe, Aerosol Sci. Tech. 4, 269)
4433    beta = 1.0_wp + knud * ( 1.142_wp + 0.558_wp * EXP( -0.999_wp / knud ) )
4434!
4435!-- 3) Particle properties
4436!-- Diffusion coefficient (m2/s) (Jacobson (2005) eq. 15.29)
4437    dfpart = beta * abo * temp / ( 3.0_wp * pi * visc * diam )
4438!
4439!-- Mean thermal velocity (m/s) (Jacobson (2005) eq. 15.32)
4440    mtvel = SQRT( ( 8.0_wp * abo * temp ) / ( pi * mpart ) )
4441!
4442!-- Particle mean free path (m) (Jacobson (2005) eq. 15.34 )
4443    omega = 8.0_wp * dfpart / ( pi * mtvel )
4444!
4445!-- Mean diameter (m)
4446    mdiam = 0.5_wp * ( diam(1) + diam(2) )
4447!
4448!-- 4) Calculation of fluxes (Brownian collision kernels) and flux matching
4449!-- following Jacobson (2005):
4450!
4451!-- Flux in continuum regime (m3/s) (eq. 15.28)
4452    flux(1) = 4.0_wp * pi * mdiam * ( dfpart(1) + dfpart(2) )
4453!
4454!-- Flux in free molec. regime (m3/s) (eq. 15.31)
4455    flux(2) = pi * SQRT( ( mtvel(1)**2 ) + ( mtvel(2)**2 ) ) * ( mdiam**2 )
4456!
4457!-- temporary variables (m) to calculate flux matching distance (m)
4458    tva(1) = ( ( mdiam + omega(1) )**3 - ( mdiam**2 + omega(1)**2 ) * SQRT( ( mdiam**2 +           &
4459               omega(1)**2 ) ) ) / ( 3.0_wp * mdiam * omega(1) ) - mdiam
4460    tva(2) = ( ( mdiam + omega(2) )**3 - ( mdiam**2 + omega(2)**2 ) * SQRT( ( mdiam**2 +           &
4461               omega(2)**2 ) ) ) / ( 3.0_wp * mdiam * omega(2) ) - mdiam
4462!
4463!-- Flux matching distance (m): the mean distance from the centre of a sphere reached by particles
4464!-- that leave sphere's surface and travel a distance of particle mean free path (eq. 15.34)
4465    fmdist = SQRT( tva(1)**2 + tva(2)**2 )
4466!
4467!-- 5) Coagulation coefficient = coalescence efficiency * collision kernel (m3/s) (eq. 15.33).
4468!--    Here assumed coalescence efficiency 1!!
4469    coagc = flux(1) / ( mdiam / ( mdiam + fmdist) + flux(1) / flux(2) )
4470!
4471!-- Corrected collision kernel (Karl et al., 2016 (ACP)): Include van der Waals and viscous forces
4472    IF ( van_der_waals_coagc )  THEN
4473       knud_p = SQRT( omega(1)**2 + omega(2)**2 ) / mdiam
4474       IF ( knud_p >= 0.1_wp  .AND.  knud_p <= 10.0_wp )  THEN
4475          coagc = coagc * ( 2.0_wp + 0.4_wp * LOG( knud_p ) )
4476       ELSE
4477          coagc = coagc * 3.0_wp
4478       ENDIF
4479    ENDIF
4480
4481 END FUNCTION coagc
4482
4483!------------------------------------------------------------------------------!
4484! Description:
4485! ------------
4486!> Calculates the change in particle volume and gas phase
4487!> concentrations due to nucleation, condensation and dissolutional growth.
4488!
4489!> Sulphuric acid and organic vapour: only condensation and no evaporation.
4490!
4491!> New gas and aerosol phase concentrations calculated according to Jacobson
4492!> (1997): Numerical techniques to solve condensational and dissolutional growth
4493!> equations when growth is coupled to reversible reactions, Aerosol Sci. Tech.,
4494!> 27, pp 491-498.
4495!
4496!> Following parameterization has been used:
4497!> Molecular diffusion coefficient of condensing vapour (m2/s)
4498!> (Reid et al. (1987): Properties of gases and liquids, McGraw-Hill, New York.)
4499!> D = {1.d-7*sqrt(1/M_air + 1/M_gas)*T^1.75} / &
4500!      {p_atm/p_stand * (d_air^(1/3) + d_gas^(1/3))^2 }
4501!> M_air = 28.965 : molar mass of air (g/mol)
4502!> d_air = 19.70  : diffusion volume of air
4503!> M_h2so4 = 98.08 : molar mass of h2so4 (g/mol)
4504!> d_h2so4 = 51.96  : diffusion volume of h2so4
4505!
4506!> Called from main aerosol model
4507!> For equations, see Jacobson, Fundamentals of Atmospheric Modeling, 2nd Edition (2005)
4508!
4509!> Coded by:
4510!> Hannele Korhonen (FMI) 2005
4511!> Harri Kokkola (FMI) 2006
4512!> Juha Tonttila (FMI) 2014
4513!> Rewritten to PALM by Mona Kurppa (UHel) 2017
4514!------------------------------------------------------------------------------!
4515 SUBROUTINE condensation( paero, pc_sa, pc_ocnv, pcocsv, pchno3, pc_nh3, pcw, pcs, ptemp, ppres,   &
4516                          ptstep, prtcl )
4517
4518    IMPLICIT NONE
4519
4520    INTEGER(iwp) ::  ss      !< start index
4521    INTEGER(iwp) ::  ee      !< end index
4522
4523    REAL(wp) ::  zcs_ocnv    !< condensation sink of nonvolatile organics (1/s)
4524    REAL(wp) ::  zcs_ocsv    !< condensation sink of semivolatile organics (1/s)
4525    REAL(wp) ::  zcs_su      !< condensation sink of sulfate (1/s)
4526    REAL(wp) ::  zcs_tot     !< total condensation sink (1/s) (gases)
4527    REAL(wp) ::  zcvap_new1  !< vapour concentration after time step (#/m3): sulphuric acid
4528    REAL(wp) ::  zcvap_new2  !< nonvolatile organics
4529    REAL(wp) ::  zcvap_new3  !< semivolatile organics
4530    REAL(wp) ::  zdfvap      !< air diffusion coefficient (m2/s)
4531    REAL(wp) ::  zdvap1      !< change in vapour concentration (#/m3): sulphuric acid
4532    REAL(wp) ::  zdvap2      !< nonvolatile organics
4533    REAL(wp) ::  zdvap3      !< semivolatile organics
4534    REAL(wp) ::  zmfp        !< mean free path of condensing vapour (m)
4535    REAL(wp) ::  zrh         !< Relative humidity [0-1]
4536    REAL(wp) ::  zvisc       !< viscosity of air (kg/(m s))
4537    REAL(wp) ::  zn_vs_c     !< ratio of nucleation of all mass transfer in the smallest bin
4538    REAL(wp) ::  zxocnv      !< ratio of organic vapour in 3nm particles
4539    REAL(wp) ::  zxsa        !< Ratio in 3nm particles: sulphuric acid
4540
4541    REAL(wp), INTENT(in) ::  ppres   !< ambient pressure (Pa)
4542    REAL(wp), INTENT(in) ::  pcs     !< Water vapour saturation concentration (kg/m3)
4543    REAL(wp), INTENT(in) ::  ptemp   !< ambient temperature (K)
4544    REAL(wp), INTENT(in) ::  ptstep  !< timestep (s)
4545
4546    REAL(wp), INTENT(inout) ::  pchno3   !< Gas concentrations (#/m3): nitric acid HNO3
4547    REAL(wp), INTENT(inout) ::  pc_nh3   !< ammonia NH3
4548    REAL(wp), INTENT(inout) ::  pc_ocnv  !< non-volatile organics
4549    REAL(wp), INTENT(inout) ::  pcocsv   !< semi-volatile organics
4550    REAL(wp), INTENT(inout) ::  pc_sa    !< sulphuric acid H2SO4
4551    REAL(wp), INTENT(inout) ::  pcw      !< Water vapor concentration (kg/m3)
4552
4553    REAL(wp), DIMENSION(nbins_aerosol)       ::  zbeta          !< transitional correction factor
4554    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate       !< collision rate (1/s)
4555    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate_ocnv  !< collision rate of OCNV (1/s)
4556    REAL(wp), DIMENSION(start_subrange_1a+1) ::  zdfpart        !< particle diffusion coef. (m2/s)
4557    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvoloc        !< change of organics volume
4558    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvolsa        !< change of sulphate volume
4559    REAL(wp), DIMENSION(2)                   ::  zj3n3          !< Formation massrate of molecules
4560                                                                !< in nucleation, (molec/m3s),
4561                                                                !< 1: H2SO4 and 2: organic vapor
4562    REAL(wp), DIMENSION(nbins_aerosol)       ::  zknud          !< particle Knudsen number
4563
4564    TYPE(component_index), INTENT(in) :: prtcl  !< Keeps track which substances are used
4565
4566    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< Aerosol properties
4567
4568    zj3n3  = 0.0_wp
4569    zrh    = pcw / pcs
4570    zxocnv = 0.0_wp
4571    zxsa   = 0.0_wp
4572!
4573!-- Nucleation
4574    IF ( nsnucl > 0 )  THEN
4575       CALL nucleation( paero, ptemp, zrh, ppres, pc_sa, pc_ocnv, pc_nh3, ptstep, zj3n3, zxsa,     &
4576                        zxocnv )
4577    ENDIF
4578!
4579!-- Condensation on pre-existing particles
4580    IF ( lscndgas )  THEN
4581!
4582!--    Initialise:
4583       zdvolsa = 0.0_wp
4584       zdvoloc = 0.0_wp
4585       zcolrate = 0.0_wp
4586!
4587!--    1) Properties of air and condensing gases:
4588!--    Viscosity of air (kg/(m s)) (Eq. 4.54 in Jabonson (2005))
4589       zvisc = ( 7.44523E-3_wp * ptemp ** 1.5_wp ) / ( 5093.0_wp * ( ptemp + 110.4_wp ) )
4590!
4591!--    Diffusion coefficient of air (m2/s)
4592       zdfvap = 5.1111E-10_wp * ptemp ** 1.75_wp * ( p_0 + 1325.0_wp ) / ppres
4593!
4594!--    Mean free path (m): same for H2SO4 and organic compounds
4595       zmfp = 3.0_wp * zdfvap * SQRT( pi * amh2so4 / ( 8.0_wp * argas * ptemp ) )
4596!
4597!--    2) Transition regime correction factor zbeta for particles (Fuchs and Sutugin (1971)):
4598!--       Size of condensing molecule considered only for nucleation mode (3 - 20 nm).
4599!
4600!--    Particle Knudsen number: condensation of gases on aerosols
4601       ss = start_subrange_1a
4602       ee = start_subrange_1a+1
4603       zknud(ss:ee) = 2.0_wp * zmfp / ( paero(ss:ee)%dwet + d_sa )
4604       ss = start_subrange_1a+2
4605       ee = end_subrange_2b
4606       zknud(ss:ee) = 2.0_wp * zmfp / paero(ss:ee)%dwet
4607!
4608!--    Transitional correction factor: aerosol + gas (the semi-empirical Fuchs- Sutugin
4609!--    interpolation function (Fuchs and Sutugin, 1971))
4610       zbeta = ( zknud + 1.0_wp ) / ( 0.377_wp * zknud + 1.0_wp + 4.0_wp / ( 3.0_wp * massacc ) *  &
4611               ( zknud + zknud ** 2 ) )
4612!
4613!--    3) Collision rate of molecules to particles
4614!--       Particle diffusion coefficient considered only for nucleation mode (3 - 20 nm)
4615!
4616!--    Particle diffusion coefficient (m2/s) (e.g. Eq. 15.29 in Jacobson (2005))
4617       zdfpart = abo * ptemp * zbeta(start_subrange_1a:start_subrange_1a+1) / ( 3.0_wp * pi * zvisc&
4618                 * paero(start_subrange_1a:start_subrange_1a+1)%dwet)
4619!
4620!--    Collision rate (mass-transfer coefficient): gases on aerosols (1/s) (Eq. 16.64 in
4621!--    Jacobson (2005))
4622       ss = start_subrange_1a
4623       ee = start_subrange_1a+1
4624       zcolrate(ss:ee) = MERGE( 2.0_wp * pi * ( paero(ss:ee)%dwet + d_sa ) * ( zdfvap + zdfpart ) *&
4625                               zbeta(ss:ee) * paero(ss:ee)%numc, 0.0_wp, paero(ss:ee)%numc > nclim )
4626       ss = start_subrange_1a+2
4627       ee = end_subrange_2b
4628       zcolrate(ss:ee) = MERGE( 2.0_wp * pi * paero(ss:ee)%dwet * zdfvap * zbeta(ss:ee) *          &
4629                                paero(ss:ee)%numc, 0.0_wp, paero(ss:ee)%numc > nclim )
4630!
4631!-- 4) Condensation sink (1/s)
4632       zcs_tot = SUM( zcolrate )   ! total sink
4633!
4634!--    5) Changes in gas-phase concentrations and particle volume
4635!
4636!--    5.1) Organic vapours
4637!
4638!--    5.1.1) Non-volatile organic compound: condenses onto all bins
4639       IF ( pc_ocnv > 1.0E+10_wp  .AND.  zcs_tot > 1.0E-30_wp  .AND. index_oc > 0 )  &
4640       THEN
4641!--       Ratio of nucleation vs. condensation rates in the smallest bin
4642          zn_vs_c = 0.0_wp
4643          IF ( zj3n3(2) > 1.0_wp )  THEN
4644             zn_vs_c = ( zj3n3(2) ) / ( zj3n3(2) + pc_ocnv * zcolrate(start_subrange_1a) )
4645          ENDIF
4646!
4647!--       Collision rate in the smallest bin, including nucleation and condensation (see
4648!--       Jacobson (2005), eq. (16.73) )
4649          zcolrate_ocnv = zcolrate
4650          zcolrate_ocnv(start_subrange_1a) = zcolrate_ocnv(start_subrange_1a) + zj3n3(2) / pc_ocnv
4651!
4652!--       Total sink for organic vapor
4653          zcs_ocnv = zcs_tot + zj3n3(2) / pc_ocnv
4654!
4655!--       New gas phase concentration (#/m3)
4656          zcvap_new2 = pc_ocnv / ( 1.0_wp + ptstep * zcs_ocnv )
4657!
4658!--       Change in gas concentration (#/m3)
4659          zdvap2 = pc_ocnv - zcvap_new2
4660!
4661!--       Updated vapour concentration (#/m3)
4662          pc_ocnv = zcvap_new2
4663!
4664!--       Volume change of particles (m3(OC)/m3(air))
4665          zdvoloc = zcolrate_ocnv(start_subrange_1a:end_subrange_2b) / zcs_ocnv * amvoc * zdvap2
4666!
4667!--       Change of volume due to condensation in 1a-2b
4668          paero(start_subrange_1a:end_subrange_2b)%volc(2) =                                       &
4669                                          paero(start_subrange_1a:end_subrange_2b)%volc(2) + zdvoloc
4670!
4671!--       Change of number concentration in the smallest bin caused by nucleation (Jacobson (2005),
4672!--       eq. (16.75)). If zxocnv = 0, then the chosen nucleation mechanism doesn't take into
4673!--       account the non-volatile organic vapors and thus the paero doesn't have to be updated.
4674          IF ( zxocnv > 0.0_wp )  THEN
4675             paero(start_subrange_1a)%numc = paero(start_subrange_1a)%numc + zn_vs_c *             &
4676                                             zdvoloc(start_subrange_1a) / amvoc / ( n3 * zxocnv )
4677          ENDIF
4678       ENDIF
4679!
4680!--    5.1.2) Semivolatile organic compound: all bins except subrange 1
4681       zcs_ocsv = SUM( zcolrate(start_subrange_2a:end_subrange_2b) ) !< sink for semi-volatile org.
4682       IF ( pcocsv > 1.0E+10_wp  .AND.  zcs_ocsv > 1.0E-30  .AND. is_used( prtcl,'OC') )  THEN
4683!
4684!--       New gas phase concentration (#/m3)
4685          zcvap_new3 = pcocsv / ( 1.0_wp + ptstep * zcs_ocsv )
4686!
4687!--       Change in gas concentration (#/m3)
4688          zdvap3 = pcocsv - zcvap_new3
4689!
4690!--       Updated gas concentration (#/m3)
4691          pcocsv = zcvap_new3
4692!
4693!--       Volume change of particles (m3(OC)/m3(air))
4694          ss = start_subrange_2a
4695          ee = end_subrange_2b
4696          zdvoloc(ss:ee) = zdvoloc(ss:ee) + zcolrate(ss:ee) / zcs_ocsv * amvoc * zdvap3
4697!
4698!--       Change of volume due to condensation in 1a-2b
4699          paero(start_subrange_1a:end_subrange_2b)%volc(2) =                                       &
4700                                          paero(start_subrange_1a:end_subrange_2b)%volc(2) + zdvoloc
4701       ENDIF
4702!
4703!--    5.2) Sulphate: condensed on all bins
4704       IF ( pc_sa > 1.0E+10_wp  .AND.  zcs_tot > 1.0E-30_wp  .AND.  index_so4 > 0 )  THEN
4705!
4706!--    Ratio of mass transfer between nucleation and condensation
4707          zn_vs_c = 0.0_wp
4708          IF ( zj3n3(1) > 1.0_wp )  THEN
4709             zn_vs_c = ( zj3n3(1) ) / ( zj3n3(1) + pc_sa * zcolrate(start_subrange_1a) )
4710          ENDIF
4711!
4712!--       Collision rate in the smallest bin, including nucleation and condensation (see
4713!--       Jacobson (2005), eq. (16.73))
4714          zcolrate(start_subrange_1a) = zcolrate(start_subrange_1a) + zj3n3(1) / pc_sa
4715!
4716!--       Total sink for sulfate (1/s)
4717          zcs_su = zcs_tot + zj3n3(1) / pc_sa
4718!
4719!--       Sulphuric acid:
4720!--       New gas phase concentration (#/m3)
4721          zcvap_new1 = pc_sa / ( 1.0_wp + ptstep * zcs_su )
4722!
4723!--       Change in gas concentration (#/m3)
4724          zdvap1 = pc_sa - zcvap_new1
4725!
4726!--       Updating vapour concentration (#/m3)
4727          pc_sa = zcvap_new1
4728!
4729!--       Volume change of particles (m3(SO4)/m3(air)) by condensation
4730          zdvolsa = zcolrate(start_subrange_1a:end_subrange_2b) / zcs_su * amvh2so4 * zdvap1
4731!
4732!--       Change of volume concentration of sulphate in aerosol [fxm]
4733          paero(start_subrange_1a:end_subrange_2b)%volc(1) =                                       &
4734                                          paero(start_subrange_1a:end_subrange_2b)%volc(1) + zdvolsa
4735!
4736!--       Change of number concentration in the smallest bin caused by nucleation
4737!--       (Jacobson (2005), equation (16.75))
4738