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

Last change on this file since 4226 was 4226, checked in by suehring, 21 months ago

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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