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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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