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

Last change on this file since 4268 was 4268, checked in by schwenkel, 20 months ago

Introducing module interface for boundary conditions and move module specific boundary conditions into their modules

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