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

Last change on this file since 4329 was 4329, checked in by motisi, 16 months ago

Renamed wall_flags_0 to wall_flags_static_0

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