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

Last change on this file since 4273 was 4273, checked in by monakurppa, 20 months ago

Add logical switched nesting_chem and nesting_offline_chem

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