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

Last change on this file since 4442 was 4442, checked in by suehring, 4 years ago

last commit documented

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