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

Last change on this file since 4118 was 4118, checked in by suehring, 2 years ago

Initialization of soil temperature and moisture via dynamic input file only for vegetation and pavement surfaces

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