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

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

last changes documented

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