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

Last change on this file since 4512 was 4512, checked in by monakurppa, 5 years ago

bug fix in salsa_mod: segmentation fault in initialising chemical components when all possible chemical components were used

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