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

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

Bug corrected in salsa_mod, salsa test case introduced and test_salsa input files updated.

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