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

Last change on this file since 4272 was 4272, checked in by schwenkel, 20 months ago

further modularization of boundary conditions: moving boundary conditions to their modules

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