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

Last change on this file since 4117 was 4117, checked in by monakurppa, 2 years ago

Apply decycling flags for salsa and set decycling boundary conditions only at the ghost points not at the prognostic grid points. See changeset 4109 for more details.

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