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

Last change on this file since 4315 was 4315, checked in by monakurppa, 22 months ago

salsa_mod: add additional check for time in PIDS_SALSA and update the urban_environment_salsa_salsa input file in the test case urban_environment_salsa

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