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

Last change on this file since 4281 was 4280, checked in by monakurppa, 19 months ago

Correct recent bugs in salsa_mod and remove chemistry specific stuff in netcdf_data_iput_mod

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