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

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

Remove salsa calls from prognostic_equations and correct a bug in the salsa deposition for urban and land surface models

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