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

Last change on this file since 4364 was 4364, checked in by monakurppa, 15 months ago

Time in the input data set relative to the start of the simulation

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