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

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

document previous changes

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