source: palm/trunk/SOURCE/salsa_mod.f90

Last change on this file was 4895, checked in by suehring, 3 months ago

Remove offset in terrain-following masked output and allow only mask_k_over_surface >= 1

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