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

Last change on this file since 4178 was 4167, checked in by suehring, 2 years ago

Merge from branch resler: Changed behaviour of masked output over surface to follow terrain and ignore buildings

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