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

Last change on this file since 4256 was 4256, checked in by monakurppa, 21 months ago

update test case urban_environment_salsa

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