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

Last change on this file since 4380 was 4380, checked in by monakurppa, 21 months ago

salsa_mod: minor bug in dry deposition and formatting

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