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

Last change on this file since 4346 was 4346, checked in by motisi, 22 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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