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

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

major changes in salsa: data input, format and performance

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