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

Last change on this file since 4313 was 4298, checked in by suehring, 19 months ago

Bugfixes: Calculation of 2-m temperature adjusted to the case the 2-m level is above the first grid level; salsa: close netcdf input files after reading; open netcdf input files with read-only attribute instead of write attribute

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