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

Last change on this file since 4342 was 4342, checked in by Giersch, 16 months ago

last commit corrected by replacing read_from_file_3d with read_from_file in all relevant p3d files, minor changes in plant_canopy_model_mod (use statements moved to module level, ocean dependency removed, redundant variables removed)

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