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

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

Initialization of soil temperature and moisture via dynamic input file only for vegetation and pavement surfaces

  • Property svn:keywords set to Id
File size: 511.2 KB
RevLine 
[2505]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!
[3885]17! Copyright 2018-2019 University of Helsinki
[3655]18! Copyright 1997-2019 Leibniz Universitaet Hannover
[2505]19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
[4118]23!
24!
25! Former revisions:
26! -----------------
27! $Id: salsa_mod.f90 4118 2019-07-25 16:11:45Z suehring $
[4117]28! - When Dirichlet condition is applied in decycling, the boundary conditions are
29!   only set at the ghost points and not at the prognostic grid points as done
30!   before
31! - Rename decycle_ns/lr to decycle_salsa_ns/lr and decycle_method to
32!   decycle_method_salsa
33! - Allocation and initialization of special advection flags salsa_advc_flags_s
34!   used for salsa. These are exclusively used for salsa variables to
35!   distinguish from the usually-used flags which might be different when
36!   decycling is applied in combination with cyclic boundary conditions.
37!   Moreover, salsa_advc_flags_s considers extended zones around buildings where
38!   the first-order upwind scheme is applied for the horizontal advection terms.
[4118]39!   This is done to overcome high concentration peaks due to stationary numerical
[4117]40!   oscillations caused by horizontal advection discretization.
[3589]41!
[4118]42! 4117 2019-07-25 08:54:02Z monakurppa
[4110]43! Pass integer flag array as well as boundary flags to WS scalar advection
44! routine
45!
46! 4109 2019-07-22 17:00:34Z suehring
[4102]47! Slightly revise setting of boundary conditions at horizontal walls, use
48! data-structure offset index instead of pre-calculate it for each facing
49!
50! 4079 2019-07-09 18:04:41Z suehring
[4079]51! Application of monotonic flux limiter for the vertical scalar advection
52! up to the topography top (only for the cache-optimized version at the
53! moment).
54!
55! 4069 2019-07-01 14:05:51Z Giersch
[4069]56! Masked output running index mid has been introduced as a local variable to
57! avoid runtime error (Loop variable has been modified) in time_integration
58!
59! 4058 2019-06-27 15:25:42Z knoop
[4058]60! Bugfix: to_be_resorted was uninitialized in case of s_H2O in 3d_data_averaging
61!
62! 4012 2019-05-31 15:19:05Z monakurppa
[4012]63! Merge salsa branch to trunk. List of changes:
64! - Error corrected in distr_update that resulted in the aerosol number size
65!   distribution not converging if the concentration was nclim.
66! - Added a separate output for aerosol liquid water (s_H2O)
67! - aerosol processes for a size bin are now calculated only if the aerosol
68!   number of concentration of that bin is > 2*nclim
69! - An initialisation error in the subroutine "deposition" corrected and the
70!   subroutine reformatted.
71! - stuff from salsa_util_mod.f90 moved into salsa_mod.f90
72! - calls for closing the netcdf input files added
73!
74! 3956 2019-05-07 12:32:52Z monakurppa
[3956]75! - Conceptual bug in depo_surf correct for urban and land surface model
76! - Subroutine salsa_tendency_ij optimized.
77! - Interfaces salsa_non_advective_processes and salsa_exchange_horiz_bounds
78!   created. These are now called in module_interface.
79!   salsa_exchange_horiz_bounds after calling salsa_driver only when needed
80!   (i.e. every dt_salsa).
81!
82! 3924 2019-04-23 09:33:06Z monakurppa
[3924]83! Correct a bug introduced by the previous update.
84!
85! 3899 2019-04-16 14:05:27Z monakurppa
[3899]86! - remove unnecessary error / location messages
87! - corrected some error message numbers
88! - allocate source arrays only if emissions or dry deposition is applied.
[3924]89!
[3899]90! 3885 2019-04-11 11:29:34Z kanani
[3885]91! Changes related to global restructuring of location messages and introduction
92! of additional debug messages
93!
94! 3876 2019-04-08 18:41:49Z knoop
[3872]95! Introduced salsa_actions module interface
96!
97! 3871 2019-04-08 14:38:39Z knoop
[3864]98! Major changes in formatting, performance and data input structure (see branch
99! the history for details)
100! - Time-dependent emissions enabled: lod=1 for yearly PM emissions that are
101!   normalised depending on the time, and lod=2 for preprocessed emissions
102!   (similar to the chemistry module).
103! - Additionally, 'uniform' emissions allowed. This emission is set constant on
104!   all horisontal upward facing surfaces and it is created based on parameters
105!   surface_aerosol_flux, aerosol_flux_dpg/sigmag/mass_fracs_a/mass_fracs_b.
106! - All emissions are now implemented as surface fluxes! No 3D sources anymore.
107! - Update the emission information by calling salsa_emission_update if
108!   skip_time_do_salsa >= time_since_reference_point and
109!   next_aero_emission_update <= time_since_reference_point
110! - Aerosol background concentrations read from PIDS_DYNAMIC. The vertical grid
111!   must match the one applied in the model.
112! - Gas emissions and background concentrations can be also read in in salsa_mod
113!   if the chemistry module is not applied.
114! - In deposition, information on the land use type can be now imported from
115!   the land use model
116! - Use SI units in PARIN, i.e. n_lognorm given in #/m3 and dpg in metres.
117! - Apply 100 character line limit
118! - Change all variable names from capital to lowercase letter
119! - Change real exponents to integer if possible. If not, precalculate the value
120!   value of exponent
121! - Rename in1a to start_subrange_1a, fn2a to end_subrange_1a etc.
122! - Rename nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and ngast -->
123!   ngases_salsa
124! - Rename ibc to index_bc, idu to index_du etc.
125! - Renamed loop indices b, c and sg to ib, ic and ig
126! - run_salsa subroutine removed
127! - Corrected a bud in salsa_driver: falsely applied ino instead of inh
128! - Call salsa_tendency within salsa_prognostic_equations which is called in
129!   module_interface_mod instead of prognostic_equations_mod
130! - Removed tailing white spaces and unused variables
131! - Change error message to start by PA instead of SA
132!
133! 3833 2019-03-28 15:04:04Z forkel
[3833]134! added USE chem_gasphase_mod for nvar, nspec and spc_names
135!
136! 3787 2019-03-07 08:43:54Z raasch
[3787]137! unused variables removed
138!
139! 3780 2019-03-05 11:19:45Z forkel
[3767]140! unused variable for file index removed from rrd-subroutines parameter list
141!
142! 3685 2019-01-21 01:02:11Z knoop
[3685]143! Some interface calls moved to module_interface + cleanup
144!
145! 3655 2019-01-07 16:51:22Z knoop
[3637]146! Implementation of the PALM module interface
147!
148! 3636 2018-12-19 13:48:34Z raasch
[3636]149! nopointer option removed
150!
151! 3630 2018-12-17 11:04:17Z knoop
[3582]152! - Moved the control parameter "salsa" from salsa_mod.f90 to control_parameters
153! - Updated salsa_rrd_local and salsa_wrd_local
154! - Add target attribute
155! - Revise initialization in case of restarts
156! - Revise masked data output
157!
[3589]158! 3582 2018-11-29 19:16:36Z suehring
[3524]159! missing comma separator inserted
160!
161! 3483 2018-11-02 14:19:26Z raasch
[3483]162! bugfix: directives added to allow compilation without netCDF
163!
164! 3481 2018-11-02 09:14:13Z raasch
[3481]165! temporary variable cc introduced to circumvent a possible Intel18 compiler bug
166! related to contiguous/non-contguous pointer/target attributes
167!
168! 3473 2018-10-30 20:50:15Z suehring
[3473]169! NetCDF input routine renamed
170!
171! 3467 2018-10-30 19:05:21Z suehring
[2525]172! Initial revision
173!
[3467]174! 3412 2018-10-24 07:25:57Z monakurppa
[2505]175!
176! Authors:
177! --------
[3582]178! @author Mona Kurppa (University of Helsinki)
[2505]179!
[3864]180!
[2505]181! Description:
182! ------------
[3864]183!> Sectional aerosol module for large scale applications SALSA
[2505]184!> (Kokkola et al., 2008, ACP 8, 2469-2483). Solves the aerosol number and mass
[3864]185!> concentration as well as chemical composition. Includes aerosol dynamic
[2505]186!> processes: nucleation, condensation/evaporation of vapours, coagulation and
[3864]187!> deposition on tree leaves, ground and roofs.
[2505]188!> Implementation is based on formulations implemented in UCLALES-SALSA except
[3864]189!> for deposition which is based on parametrisations by Zhang et al. (2001,
190!> Atmos. Environ. 35, 549-560) or Petroff&Zhang (2010, Geosci. Model Dev. 3,
[2505]191!> 753-769)
192!>
[3864]193!> @todo Apply information from emission_stack_height to lift emission sources
194!> @todo emission mode "parameterized", i.e. based on street type
[3899]195!> @todo Allow insoluble emissions
[4012]196!> @todo Apply flux limiter in prognostic equations
[2505]197!------------------------------------------------------------------------------!
198 MODULE salsa_mod
[3308]199
[4012]200    USE basic_constants_and_equations_mod,                                                         &
[3308]201        ONLY:  c_p, g, p_0, pi, r_d
[3864]202
[4012]203    USE chem_gasphase_mod,                                                                         &
[3833]204        ONLY:  nspec, nvar, spc_names
205
[4012]206    USE chem_modules,                                                                              &
[3876]207        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, chem_species
[2754]208
[4117]209    USE control_parameters,                                                                        &
210        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s,      &
211               bc_lr, bc_lr_cyc, bc_ns, bc_ns_cyc, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
212               bc_radiation_s, coupling_char, debug_output, dt_3d, intermediate_timestep_count,    &
213               intermediate_timestep_count_max, land_surface, message_string, monotonic_limiter_z, &
214               plant_canopy, pt_surface, salsa, scalar_advec, surface_pressure,                    &
215               time_since_reference_point, timestep_scheme, tsc, urban_surface, ws_scheme_sca
[3637]216
[4012]217    USE indices,                                                                                   &
218        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nz, nzt, wall_flags_0
[3864]219
[2505]220    USE kinds
[3864]221
[4012]222    USE netcdf_data_input_mod,                                                                     &
223        ONLY:  chem_emis_att_type, chem_emis_val_type
224
[2505]225    USE pegrid
[3864]226
[4012]227    USE statistics,                                                                                &
[3871]228        ONLY:  sums_salsa_ws_l
229
[2505]230    IMPLICIT NONE
231!
232!-- SALSA constants:
233!
234!-- Local constants:
[3899]235    INTEGER(iwp), PARAMETER ::  luc_urban = 15     !< default landuse type for urban
[3956]236    INTEGER(iwp), PARAMETER ::  ngases_salsa  = 5  !< total number of gaseous tracers:
[3864]237                                                   !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV
238                                                   !< (non-volatile OC), 5 = OCSV (semi-volatile)
[3956]239    INTEGER(iwp), PARAMETER ::  nmod = 7     !< number of modes for initialising the aerosol size
[3864]240                                             !< distribution
[3956]241    INTEGER(iwp), PARAMETER ::  nreg = 2     !< Number of main size subranges
[3864]242    INTEGER(iwp), PARAMETER ::  maxspec = 7  !< Max. number of aerosol species
243    INTEGER(iwp), PARAMETER ::  season = 1   !< For dry depostion by Zhang et al.: 1 = summer,
244                                             !< 2 = autumn (no harvest yet), 3 = late autumn
245                                             !< (already frost), 4 = winter, 5 = transitional spring
246!
[2505]247!-- Universal constants
[3864]248    REAL(wp), PARAMETER ::  abo    = 1.380662E-23_wp   !< Boltzmann constant (J/K)
249    REAL(wp), PARAMETER ::  alv    = 2.260E+6_wp       !< latent heat for H2O
250                                                       !< vaporisation (J/kg)
251    REAL(wp), PARAMETER ::  alv_d_rv  = 4896.96865_wp  !< alv / rv
252    REAL(wp), PARAMETER ::  am_airmol = 4.8096E-26_wp  !< Average mass of one air
253                                                       !< molecule (Jacobson,
254                                                       !< 2005, Eq. 2.3)
255    REAL(wp), PARAMETER ::  api6   = 0.5235988_wp      !< pi / 6
256    REAL(wp), PARAMETER ::  argas  = 8.314409_wp       !< Gas constant (J/(mol K))
257    REAL(wp), PARAMETER ::  argas_d_cpd = 8.281283865E-3_wp  !< argas per cpd
258    REAL(wp), PARAMETER ::  avo    = 6.02214E+23_wp    !< Avogadro constant (1/mol)
259    REAL(wp), PARAMETER ::  d_sa   = 5.539376964394570E-10_wp  !< diameter of condensing sulphuric
260                                                               !< acid molecule (m)
261    REAL(wp), PARAMETER ::  for_ppm_to_nconc =  7.243016311E+16_wp !< ppm * avo / R (K/(Pa*m3))
[2505]262    REAL(wp), PARAMETER ::  epsoc  = 0.15_wp          !< water uptake of organic
[3864]263                                                      !< material
264    REAL(wp), PARAMETER ::  mclim  = 1.0E-23_wp       !< mass concentration min limit (kg/m3)
265    REAL(wp), PARAMETER ::  n3     = 158.79_wp        !< Number of H2SO4 molecules in 3 nm cluster
266                                                      !< if d_sa=5.54e-10m
267    REAL(wp), PARAMETER ::  nclim  = 1.0_wp           !< number concentration min limit (#/m3)
268    REAL(wp), PARAMETER ::  surfw0 = 0.073_wp         !< surface tension of water at 293 K (J/m2)
269!
[2505]270!-- Molar masses in kg/mol
[3956]271    REAL(wp), PARAMETER ::  ambc     = 12.0E-3_wp     !< black carbon (BC)
272    REAL(wp), PARAMETER ::  amdair   = 28.970E-3_wp   !< dry air
273    REAL(wp), PARAMETER ::  amdu     = 100.E-3_wp     !< mineral dust
274    REAL(wp), PARAMETER ::  amh2o    = 18.0154E-3_wp  !< H2O
275    REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp    !< H2SO4
276    REAL(wp), PARAMETER ::  amhno3   = 63.01E-3_wp    !< HNO3
277    REAL(wp), PARAMETER ::  amn2o    = 44.013E-3_wp   !< N2O
278    REAL(wp), PARAMETER ::  amnh3    = 17.031E-3_wp   !< NH3
279    REAL(wp), PARAMETER ::  amo2     = 31.9988E-3_wp  !< O2
280    REAL(wp), PARAMETER ::  amo3     = 47.998E-3_wp   !< O3
281    REAL(wp), PARAMETER ::  amoc     = 150.E-3_wp     !< organic carbon (OC)
282    REAL(wp), PARAMETER ::  amss     = 58.44E-3_wp    !< sea salt (NaCl)
[3864]283!
[2505]284!-- Densities in kg/m3
[3864]285    REAL(wp), PARAMETER ::  arhobc     = 2000.0_wp  !< black carbon
286    REAL(wp), PARAMETER ::  arhodu     = 2650.0_wp  !< mineral dust
287    REAL(wp), PARAMETER ::  arhoh2o    = 1000.0_wp  !< H2O
288    REAL(wp), PARAMETER ::  arhoh2so4  = 1830.0_wp  !< SO4
289    REAL(wp), PARAMETER ::  arhohno3   = 1479.0_wp  !< HNO3
290    REAL(wp), PARAMETER ::  arhonh3    = 1530.0_wp  !< NH3
291    REAL(wp), PARAMETER ::  arhooc     = 2000.0_wp  !< organic carbon
292    REAL(wp), PARAMETER ::  arhoss     = 2165.0_wp  !< sea salt (NaCl)
293!
[2505]294!-- Volume of molecule in m3/#
295    REAL(wp), PARAMETER ::  amvh2o   = amh2o /avo / arhoh2o      !< H2O
296    REAL(wp), PARAMETER ::  amvh2so4 = amh2so4 / avo / arhoh2so4 !< SO4
[3864]297    REAL(wp), PARAMETER ::  amvhno3  = amhno3 / avo / arhohno3   !< HNO3
298    REAL(wp), PARAMETER ::  amvnh3   = amnh3 / avo / arhonh3     !< NH3
[2505]299    REAL(wp), PARAMETER ::  amvoc    = amoc / avo / arhooc       !< OC
300    REAL(wp), PARAMETER ::  amvss    = amss / avo / arhoss       !< sea salt
301!
[3864]302!-- Constants for the dry deposition model by Petroff and Zhang (2010):
303!-- obstacle characteristic dimension "L" (cm) (plane obstacle by default) and empirical constants
304!-- C_B, C_IN, C_IM, beta_IM and C_IT for each land use category (15, as in Zhang et al. (2001))
305    REAL(wp), DIMENSION(1:15), PARAMETER :: l_p10 = &
[3956]306        (/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/)
[3864]307    REAL(wp), DIMENSION(1:15), PARAMETER :: c_b_p10 = &
[3956]308        (/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/)
[3864]309    REAL(wp), DIMENSION(1:15), PARAMETER :: c_in_p10 = &
[3956]310        (/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/)
[3864]311    REAL(wp), DIMENSION(1:15), PARAMETER :: c_im_p10 = &
[3956]312        (/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/)
[3864]313    REAL(wp), DIMENSION(1:15), PARAMETER :: beta_im_p10 = &
[3956]314        (/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/)
[3864]315    REAL(wp), DIMENSION(1:15), PARAMETER :: c_it_p10 = &
[3956]316        (/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/)
[3864]317!
318!-- Constants for the dry deposition model by Zhang et al. (2001):
319!-- empirical constants "alpha" and "gamma" and characteristic radius "A" for
320!-- each land use category (15) and season (5)
321    REAL(wp), DIMENSION(1:15), PARAMETER :: alpha_z01 = &
322        (/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/)
323    REAL(wp), DIMENSION(1:15), PARAMETER :: gamma_z01 = &
324        (/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/)
325    REAL(wp), DIMENSION(1:15,1:5), PARAMETER :: A_z01 =  RESHAPE( (/& 
326         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
327         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
328         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
329         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
330         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
331                                                           /), (/ 15, 5 /) )
332!-- Land use categories (based on Z01 but the same applies here also for P10):
333!-- 1 = evergreen needleleaf trees,
334!-- 2 = evergreen broadleaf trees,
335!-- 3 = deciduous needleleaf trees,
336!-- 4 = deciduous broadleaf trees,
337!-- 5 = mixed broadleaf and needleleaf trees (deciduous broadleaf trees for P10),
338!-- 6 = grass (short grass for P10),
339!-- 7 = crops, mixed farming,
340!-- 8 = desert,
341!-- 9 = tundra,
342!-- 10 = shrubs and interrupted woodlands (thorn shrubs for P10),
343!-- 11 = wetland with plants (long grass for P10)
344!-- 12 = ice cap and glacier,
345!-- 13 = inland water (inland lake for P10)
346!-- 14 = ocean (water for P10),
347!-- 15 = urban
348!
349!-- SALSA variables:
350    CHARACTER(LEN=20)  ::  bc_salsa_b = 'neumann'                 !< bottom boundary condition
351    CHARACTER(LEN=20)  ::  bc_salsa_t = 'neumann'                 !< top boundary condition
352    CHARACTER(LEN=20)  ::  depo_pcm_par = 'zhang2001'             !< or 'petroff2010'
353    CHARACTER(LEN=20)  ::  depo_pcm_type = 'deciduous_broadleaf'  !< leaf type
354    CHARACTER(LEN=20)  ::  depo_surf_par = 'zhang2001'            !< or 'petroff2010'
355    CHARACTER(LEN=100) ::  input_file_dynamic = 'PIDS_DYNAMIC'    !< file name for dynamic input
356    CHARACTER(LEN=100) ::  input_file_salsa   = 'PIDS_SALSA'      !< file name for emission data
357    CHARACTER(LEN=20)  ::  salsa_emission_mode = 'no_emission'    !< 'no_emission', 'uniform',
358                                                                  !< 'parameterized', 'read_from_file'
359
[4117]360    CHARACTER(LEN=20), DIMENSION(4) ::  decycle_method_salsa =                                     &
[3864]361                                                 (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
362                                     !< Decycling method at horizontal boundaries
363                                     !< 1=left, 2=right, 3=south, 4=north
364                                     !< dirichlet = initial profiles for the ghost and first 3 layers
365                                     !< neumann = zero gradient
366
367    CHARACTER(LEN=3), DIMENSION(maxspec) ::  listspec = &  !< Active aerosols
368                                   (/'SO4','   ','   ','   ','   ','   ','   '/)
369
[3956]370    INTEGER(iwp) ::  depo_pcm_par_num = 1   !< parametrisation type: 1=zhang2001, 2=petroff2010
[3864]371    INTEGER(iwp) ::  depo_pcm_type_num = 0  !< index for the dry deposition type on the plant canopy
[3956]372    INTEGER(iwp) ::  depo_surf_par_num = 1  !< parametrisation type: 1=zhang2001, 2=petroff2010
[3864]373    INTEGER(iwp) ::  dots_salsa = 0         !< starting index for salsa-timeseries
374    INTEGER(iwp) ::  end_subrange_1a = 1    !< last index for bin subrange 1a
375    INTEGER(iwp) ::  end_subrange_2a = 1    !< last index for bin subrange 2a
376    INTEGER(iwp) ::  end_subrange_2b = 1    !< last index for bin subrange 2b
377    INTEGER(iwp) ::  ibc_salsa_b            !< index for the bottom boundary condition
378    INTEGER(iwp) ::  ibc_salsa_t            !< index for the top boundary condition
379    INTEGER(iwp) ::  index_bc  = -1         !< index for black carbon (BC)
380    INTEGER(iwp) ::  index_du  = -1         !< index for dust
381    INTEGER(iwp) ::  index_nh  = -1         !< index for NH3
382    INTEGER(iwp) ::  index_no  = -1         !< index for HNO3
383    INTEGER(iwp) ::  index_oc  = -1         !< index for organic carbon (OC)
[3899]384    INTEGER(iwp) ::  index_so4 = -1         !< index for SO4 or H2SO4
385    INTEGER(iwp) ::  index_ss  = -1         !< index for sea salt
386    INTEGER(iwp) ::  init_aerosol_type = 0  !< Initial size distribution type
[3864]387                                            !< 0 = uniform (read from PARIN)
388                                            !< 1 = read vertical profile of the mode number
389                                            !<     concentration from an input file
[3899]390    INTEGER(iwp) ::  init_gases_type = 0    !< Initial gas concentration type
391                                            !< 0 = uniform (read from PARIN)
392                                            !< 1 = read vertical profile from an input file
[3864]393    INTEGER(iwp) ::  lod_gas_emissions = 0  !< level of detail of the gaseous emission data
394    INTEGER(iwp) ::  nbins_aerosol = 1      !< total number of size bins
395    INTEGER(iwp) ::  ncc   = 1              !< number of chemical components used
396    INTEGER(iwp) ::  ncomponents_mass = 1   !< total number of chemical compounds (ncc+1)
397                                            !< if particle water is advected)
398    INTEGER(iwp) ::  nj3 = 1                !< J3 parametrization (nucleation)
399                                            !< 1 = condensational sink (Kerminen&Kulmala, 2002)
400                                            !< 2 = coagulational sink (Lehtinen et al. 2007)
401                                            !< 3 = coagS+self-coagulation (Anttila et al. 2010)
402    INTEGER(iwp) ::  nsnucl = 0             !< Choice of the nucleation scheme:
403                                            !< 0 = off
404                                            !< 1 = binary nucleation
405                                            !< 2 = activation type nucleation
406                                            !< 3 = kinetic nucleation
407                                            !< 4 = ternary nucleation
408                                            !< 5 = nucleation with ORGANICs
409                                            !< 6 = activation type of nucleation with H2SO4+ORG
410                                            !< 7 = heteromolecular nucleation with H2SO4*ORG
411                                            !< 8 = homomolecular nucleation of H2SO4
412                                            !<     + heteromolecular nucleation with H2SO4*ORG
413                                            !< 9 = homomolecular nucleation of H2SO4 and ORG
414                                            !<     + heteromolecular nucleation with H2SO4*ORG
415    INTEGER(iwp) ::  start_subrange_1a = 1  !< start index for bin subranges: subrange 1a
416    INTEGER(iwp) ::  start_subrange_2a = 1  !<                                subrange 2a
417    INTEGER(iwp) ::  start_subrange_2b = 1  !<                                subrange 2b
418
419    INTEGER(iwp), DIMENSION(nreg) ::  nbin = (/ 3, 7/)  !< Number of size bins per subrange: 1 & 2
420
421    INTEGER(iwp), DIMENSION(ngases_salsa) ::  gas_index_chem = &
422                                                 (/ 1, 1, 1, 1, 1/)  !< gas indices in chemistry_model_mod
423                                                 !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV, 5 = OCSV
424    INTEGER(iwp), DIMENSION(ngases_salsa) ::  emission_index_chem  !< gas indices in the gas emission file
425
426    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  k_topo_top  !< vertical index of the topography top
[4117]427
428    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE  ::  salsa_advc_flags_s !< flags used to degrade order of advection
429                                                                        !< scheme for salsa variables near walls and
430                                                                        !< lateral boundaries
[3864]431!
[2505]432!-- SALSA switches:
[3956]433    LOGICAL ::  advect_particle_water   = .TRUE.   !< Advect water concentration of particles
[4117]434    LOGICAL ::  decycle_salsa_lr        = .FALSE.  !< Undo cyclic boundaries: left and right
435    LOGICAL ::  decycle_salsa_ns        = .FALSE.  !< Undo cyclic boundaries: north and south
[3956]436    LOGICAL ::  include_emission        = .FALSE.  !< Include or not emissions
437    LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
438    LOGICAL ::  nest_salsa              = .FALSE.  !< Apply nesting for salsa
439    LOGICAL ::  no_insoluble            = .FALSE.  !< Exclude insoluble chemical components
440    LOGICAL ::  read_restart_data_salsa = .FALSE.  !< Read restart data for salsa
441    LOGICAL ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA from
442                                                   !< the chemistry model
443    LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Enhancement of coagulation kernel by van der
[3864]444                                                   !< Waals and viscous forces
[3956]445    LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
[3864]446!
[2505]447!-- Process switches: nl* is read from the NAMELIST and is NOT changed.
[3864]448!--                   ls* is the switch used and will get the value of nl*
[2505]449!--                       except for special circumstances (spinup period etc.)
[3864]450    LOGICAL ::  nlcoag       = .FALSE.  !< Coagulation master switch
451    LOGICAL ::  lscoag       = .FALSE.  !<
452    LOGICAL ::  nlcnd        = .FALSE.  !< Condensation master switch
453    LOGICAL ::  lscnd        = .FALSE.  !<
454    LOGICAL ::  nlcndgas     = .FALSE.  !< Condensation of precursor gases
455    LOGICAL ::  lscndgas     = .FALSE.  !<
456    LOGICAL ::  nlcndh2oae   = .FALSE.  !< Condensation of H2O on aerosol
457    LOGICAL ::  lscndh2oae   = .FALSE.  !< particles (FALSE -> equilibrium calc.)
458    LOGICAL ::  nldepo       = .FALSE.  !< Deposition master switch
459    LOGICAL ::  lsdepo       = .FALSE.  !<
460    LOGICAL ::  nldepo_surf  = .FALSE.  !< Deposition on vegetation master switch
461    LOGICAL ::  lsdepo_surf  = .FALSE.  !<
462    LOGICAL ::  nldepo_pcm   = .FALSE.  !< Deposition on walls master switch
463    LOGICAL ::  lsdepo_pcm   = .FALSE.  !<
464    LOGICAL ::  nldistupdate = .TRUE.   !< Size distribution update master switch
465    LOGICAL ::  lsdistupdate = .FALSE.  !<
466    LOGICAL ::  lspartition  = .FALSE.  !< Partition of HNO3 and NH3
467
468    REAL(wp) ::  act_coeff = 1.0E-7_wp               !< Activation coefficient
469    REAL(wp) ::  dt_salsa  = 0.00001_wp              !< Time step of SALSA
470    REAL(wp) ::  h2so4_init = nclim                  !< Init value for sulphuric acid gas
471    REAL(wp) ::  hno3_init  = nclim                  !< Init value for nitric acid gas
472    REAL(wp) ::  last_salsa_time = 0.0_wp            !< previous salsa call
473    REAL(wp) ::  next_aero_emission_update = 0.0_wp  !< previous emission update
474    REAL(wp) ::  next_gas_emission_update = 0.0_wp   !< previous emission update
475    REAL(wp) ::  nf2a = 1.0_wp                       !< Number fraction allocated to 2a-bins
476    REAL(wp) ::  nh3_init  = nclim                   !< Init value for ammonia gas
477    REAL(wp) ::  ocnv_init = nclim                   !< Init value for non-volatile organic gases
478    REAL(wp) ::  ocsv_init = nclim                   !< Init value for semi-volatile organic gases
479    REAL(wp) ::  rhlim = 1.20_wp                     !< RH limit in %/100. Prevents unrealistical RH
480    REAL(wp) ::  skip_time_do_salsa = 0.0_wp         !< Starting time of SALSA (s)
[2505]481!
[3864]482!-- Initial log-normal size distribution: mode diameter (dpg, metres),
483!-- standard deviation (sigmag) and concentration (n_lognorm, #/m3)
484    REAL(wp), DIMENSION(nmod) ::  dpg   = &
[3899]485                     (/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/)
[3864]486    REAL(wp), DIMENSION(nmod) ::  sigmag  = &
487                                        (/1.8_wp, 2.16_wp, 2.21_wp, 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/)
488    REAL(wp), DIMENSION(nmod) ::  n_lognorm = &
489                             (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
490!
491!-- Initial mass fractions / chemical composition of the size distribution
492    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_a = & !< mass fractions between
493             (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for A bins
494    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_b = & !< mass fractions between
495             (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for B bins
[2777]496    REAL(wp), DIMENSION(nreg+1) ::  reglim = & !< Min&max diameters of size subranges
497                                 (/ 3.0E-9_wp, 5.0E-8_wp, 1.0E-5_wp/)
[2505]498!
[3864]499!-- Initial log-normal size distribution: mode diameter (dpg, metres), standard deviation (sigmag)
500!-- concentration (n_lognorm, #/m3) and mass fractions of all chemical components (listed in
501!-- listspec) for both a (soluble) and b (insoluble) bins.
502    REAL(wp), DIMENSION(nmod) ::  aerosol_flux_dpg   = &
[3899]503                     (/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/)
[3864]504    REAL(wp), DIMENSION(nmod) ::  aerosol_flux_sigmag  = &
505                                        (/1.8_wp, 2.16_wp, 2.21_wp, 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/)
506    REAL(wp), DIMENSION(maxspec) ::  aerosol_flux_mass_fracs_a = &
507                                                               (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
508    REAL(wp), DIMENSION(maxspec) ::  aerosol_flux_mass_fracs_b = &
509                                                               (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
[3899]510    REAL(wp), DIMENSION(nmod) ::  surface_aerosol_flux = &
[4012]511                                 (/1.0E+8_wp, 1.0E+9_wp, 1.0E+5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
[3864]512
513    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bin_low_limits     !< to deliver information about
514                                                               !< the lower diameters per bin
515    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bc_am_t_val        !< vertical gradient of: aerosol mass
516    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bc_an_t_val        !< of: aerosol number
517    REAL(wp), DIMENSION(:), ALLOCATABLE ::  bc_gt_t_val        !< salsa gases near domain top
518    REAL(wp), DIMENSION(:), ALLOCATABLE ::  gas_emission_time  !< Time array in gas emission data (s)
519    REAL(wp), DIMENSION(:), ALLOCATABLE ::  nsect              !< Background number concentrations
520    REAL(wp), DIMENSION(:), ALLOCATABLE ::  massacc            !< Mass accomodation coefficients
[2525]521!
[3864]522!-- SALSA derived datatypes:
523!
[4012]524!-- Component index
525    TYPE component_index
526       CHARACTER(len=3), ALLOCATABLE ::  comp(:)  !< Component name
527       INTEGER(iwp) ::  ncomp  !< Number of components
528       INTEGER(iwp), ALLOCATABLE ::  ind(:)  !< Component index
529    END TYPE component_index
530!
[3956]531!-- For matching LSM and USM surface types and the deposition module surface types
532    TYPE match_surface
533       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_lupg  !< index for pavement / green roofs
534       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luvw  !< index for vegetation / walls
535       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luww  !< index for water / windows
536    END TYPE match_surface
[3864]537!
538!-- Aerosol emission data attributes
539    TYPE salsa_emission_attribute_type
540
541       CHARACTER(LEN=25) ::   units
542
543       CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE ::   cat_name    !<
544       CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE ::   cc_name     !<
545       CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE ::   unit_time   !<
546       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names   !<
547
548       INTEGER(iwp) ::  lod = 0            !< level of detail
549       INTEGER(iwp) ::  nbins = 10         !< number of aerosol size bins
550       INTEGER(iwp) ::  ncat  = 0          !< number of emission categories
551       INTEGER(iwp) ::  ncc   = 7          !< number of aerosol chemical components
552       INTEGER(iwp) ::  nhoursyear = 0     !< number of hours: HOURLY mode
553       INTEGER(iwp) ::  nmonthdayhour = 0  !< number of month days and hours: MDH mode
554       INTEGER(iwp) ::  num_vars           !< number of variables
555       INTEGER(iwp) ::  nt  = 0            !< number of time steps
556       INTEGER(iwp) ::  nz  = 0            !< number of vertical levels
557       INTEGER(iwp) ::  tind               !< time index for reference time in salsa emission data
558
559       INTEGER(iwp), DIMENSION(maxspec) ::  cc_input_to_model   !<
560
561       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  cat_index  !< Index of emission categories
562       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  cc_index   !< Index of chemical components
563
564       REAL(wp) ::  conversion_factor  !< unit conversion factor for aerosol emissions
565
566       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dmid         !< mean diameters of size bins (m)
567       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rho          !< average density (kg/m3)
568       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time         !< time (s)
569       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_factor  !< emission time factor
570       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z            !< height (m)
571
572       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  etf  !< emission time factor
573       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: stack_height
574
575    END TYPE salsa_emission_attribute_type
576!
577!-- The default size distribution and mass composition per emission category:
578!-- 1 = traffic, 2 = road dust, 3 = wood combustion, 4 = other
579!-- Mass fractions: H2SO4, OC, BC, DU, SS, HNO3, NH3
580    TYPE salsa_emission_mode_type
581
582       INTEGER(iwp) ::  ndm = 3  !< number of default modes
583       INTEGER(iwp) ::  ndc = 4  !< number of default categories
584
585       CHARACTER(LEN=25), DIMENSION(1:4) ::  cat_name_table = (/'traffic exhaust', &
586                                                                'road dust      ', &
587                                                                'wood combustion', &
588                                                                'other          '/)
589
590       INTEGER(iwp), DIMENSION(1:4) ::  cat_input_to_model   !<
591
592       REAL(wp), DIMENSION(1:3) ::  dpg_table = (/ 13.5E-9_wp, 1.4E-6_wp, 5.4E-8_wp/)  !<
593       REAL(wp), DIMENSION(1:3) ::  ntot_table  !<
594       REAL(wp), DIMENSION(1:3) ::  sigmag_table = (/ 1.6_wp, 1.4_wp, 1.7_wp /)  !<
595
596       REAL(wp), DIMENSION(1:maxspec,1:4) ::  mass_frac_table = &  !<
597          RESHAPE( (/ 0.04_wp, 0.48_wp, 0.48_wp, 0.0_wp,  0.0_wp, 0.0_wp, 0.0_wp, &
598                      0.0_wp,  0.05_wp, 0.0_wp,  0.95_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
599                      0.0_wp,  0.5_wp,  0.5_wp,  0.0_wp,  0.0_wp, 0.0_wp, 0.0_wp, &
600                      0.0_wp,  0.5_wp,  0.5_wp,  0.0_wp,  0.0_wp, 0.0_wp, 0.0_wp  &
601                   /), (/maxspec,4/) )
602
603       REAL(wp), DIMENSION(1:3,1:4) ::  pm_frac_table = & !< rel. mass
604                                     RESHAPE( (/ 0.016_wp, 0.000_wp, 0.984_wp, &
605                                                 0.000_wp, 1.000_wp, 0.000_wp, &
606                                                 0.000_wp, 0.000_wp, 1.000_wp, &
607                                                 1.000_wp, 0.000_wp, 1.000_wp  &
608                                              /), (/3,4/) )
609
610    END TYPE salsa_emission_mode_type
611!
612!-- Aerosol emission data values
613    TYPE salsa_emission_value_type
614
615       REAL(wp) ::  fill  !< fill value
616
617       REAL(wp), DIMENSION(:), ALLOCATABLE :: preproc_mass_fracs  !< mass fractions
618
619       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: def_mass_fracs  !< mass fractions per emis. category
620
621       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: def_data      !< surface emission values in PM
622       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: preproc_data  !< surface emission values per bin
623
624    END TYPE salsa_emission_value_type
625!
626!-- Prognostic variable: Aerosol size bin information (number (#/m3) and mass (kg/m3) concentration)
627!-- and the concentration of gaseous tracers (#/m3). Gas tracers are contained sequentially in
628!-- dimension 4 as:
629!-- 1. H2SO4, 2. HNO3, 3. NH3, 4. OCNV (non-volatile organics), 5. OCSV (semi-volatile)
[2525]630    TYPE salsa_variable
[3864]631
632       REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  init  !<
633
634       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  diss_s     !<
635       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  flux_s     !<
636       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  source     !<
637       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  sums_ws_l  !<
638
639       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l  !<
640       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l  !<
641
642       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  conc     !<
643       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  conc_p   !<
644       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  tconc_m  !<
645
[2525]646    END TYPE salsa_variable
[3864]647!
648!-- Datatype used to store information about the binned size distributions of aerosols
[2505]649    TYPE t_section
[3864]650
651       REAL(wp) ::  dmid     !< bin middle diameter (m)
[2505]652       REAL(wp) ::  vhilim   !< bin volume at the high limit
653       REAL(wp) ::  vlolim   !< bin volume at the low limit
654       REAL(wp) ::  vratiohi !< volume ratio between the center and high limit
655       REAL(wp) ::  vratiolo !< volume ratio between the center and low limit
656       !******************************************************
657       ! ^ Do NOT change the stuff above after initialization !
658       !******************************************************
[3864]659       REAL(wp) ::  core    !< Volume of dry particle
[2505]660       REAL(wp) ::  dwet    !< Wet diameter or mean droplet diameter (m)
[3864]661       REAL(wp) ::  numc    !< Number concentration of particles/droplets (#/m3)
[2505]662       REAL(wp) ::  veqh2o  !< Equilibrium H2O concentration for each particle
[3864]663
664       REAL(wp), DIMENSION(maxspec+1) ::  volc !< Volume concentrations (m^3/m^3) of aerosols +
665                                               !< water. Since most of the stuff in SALSA is hard
666                                               !< coded, these *have to be* in the order
667                                               !< 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
668    END TYPE t_section
669
670    TYPE(salsa_emission_attribute_type) ::  aero_emission_att  !< emission attributes
671    TYPE(salsa_emission_value_type)     ::  aero_emission      !< emission values
672    TYPE(salsa_emission_mode_type)      ::  def_modes          !< default emission modes
673
[4012]674    TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
675    TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
676
[3864]677    TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
678
[3956]679    TYPE(match_surface) ::  lsm_to_depo_h  !< to match the deposition module and horizontal LSM surfaces
680    TYPE(match_surface) ::  usm_to_depo_h  !< to match the deposition module and horizontal USM surfaces
[3864]681
[3956]682    TYPE(match_surface), DIMENSION(0:3) ::  lsm_to_depo_v  !< to match the deposition mod. and vertical LSM surfaces
683    TYPE(match_surface), DIMENSION(0:3) ::  usm_to_depo_v  !< to match the deposition mod. and vertical USM surfaces
[2505]684!
[3864]685!-- SALSA variables: as x = x(k,j,i,bin).
686!-- The 4th dimension contains all the size bins sequentially for each aerosol species  + water.
[2505]687!
[3864]688!-- Prognostic variables:
[2505]689!
[3864]690!-- Number concentration (#/m3)
691    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_number  !<
692    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_1  !<
693    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_2  !<
694    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_3  !<
[2505]695!
[3864]696!-- Mass concentration (kg/m3)
697    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_mass  !<
698    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_1  !<
699    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_2  !<
700    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_3  !<
[2505]701!
[3864]702!-- Gaseous concentrations (#/m3)
703    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  salsa_gas  !<
704    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_1  !<
705    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_2  !<
706    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_3  !<
[2505]707!
708!-- Diagnostic tracers
[3864]709    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  sedim_vd  !< sedimentation velocity per bin (m/s)
710    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  ra_dry    !< aerosol dry radius (m)
711
[2505]712!-- Particle component index tables
[3864]713    TYPE(component_index) :: prtcl  !< Contains "getIndex" which gives the index for a given aerosol
714                                    !< component name: 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
715!
[2505]716!-- Data output arrays:
[3864]717!
[3412]718!-- Gases:
[3864]719    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_h2so4_av  !< H2SO4
720    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_hno3_av   !< HNO3
721    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_nh3_av    !< NH3
722    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_ocnv_av   !< non-volatile OC
723    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_ocsv_av   !< semi-volatile OC
724!
725!-- Integrated:
726    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ldsa_av  !< lung-deposited surface area
727    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ntot_av  !< total number concentration
728    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  pm25_av  !< PM2.5
729    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  pm10_av  !< PM10
730!
731!-- In the particle phase:
732    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_bc_av   !< black carbon
733    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_du_av   !< dust
734    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_h2o_av  !< liquid water
735    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_nh_av   !< ammonia
736    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_no_av   !< nitrates
737    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_oc_av   !< org. carbon
738    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_so4_av  !< sulphates
739    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_ss_av   !< sea salt
740!
741!-- Bin specific mass and number concentrations:
742    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mbins_av  !< bin mas
743    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nbins_av  !< bin number
[3582]744
[2505]745!
746!-- PALM interfaces:
747!
748!-- Boundary conditions:
749    INTERFACE salsa_boundary_conds
750       MODULE PROCEDURE salsa_boundary_conds
751       MODULE PROCEDURE salsa_boundary_conds_decycle
752    END INTERFACE salsa_boundary_conds
[3864]753!
[2505]754!-- Data output checks for 2D/3D data to be done in check_parameters
755    INTERFACE salsa_check_data_output
756       MODULE PROCEDURE salsa_check_data_output
757    END INTERFACE salsa_check_data_output
758!
759!-- Input parameter checks to be done in check_parameters
760    INTERFACE salsa_check_parameters
761       MODULE PROCEDURE salsa_check_parameters
762    END INTERFACE salsa_check_parameters
763!
764!-- Averaging of 3D data for output
765    INTERFACE salsa_3d_data_averaging
766       MODULE PROCEDURE salsa_3d_data_averaging
767    END INTERFACE salsa_3d_data_averaging
768!
769!-- Data output of 2D quantities
770    INTERFACE salsa_data_output_2d
771       MODULE PROCEDURE salsa_data_output_2d
772    END INTERFACE salsa_data_output_2d
773!
774!-- Data output of 3D data
775    INTERFACE salsa_data_output_3d
776       MODULE PROCEDURE salsa_data_output_3d
777    END INTERFACE salsa_data_output_3d
778!
779!-- Data output of 3D data
780    INTERFACE salsa_data_output_mask
781       MODULE PROCEDURE salsa_data_output_mask
782    END INTERFACE salsa_data_output_mask
783!
784!-- Definition of data output quantities
785    INTERFACE salsa_define_netcdf_grid
786       MODULE PROCEDURE salsa_define_netcdf_grid
[2730]787    END INTERFACE salsa_define_netcdf_grid
[2505]788!
789!-- Output of information to the header file
790    INTERFACE salsa_header
791       MODULE PROCEDURE salsa_header
792    END INTERFACE salsa_header
793!
[3864]794!-- Initialization actions
[2505]795    INTERFACE salsa_init
796       MODULE PROCEDURE salsa_init
797    END INTERFACE salsa_init
798!
799!-- Initialization of arrays
800    INTERFACE salsa_init_arrays
801       MODULE PROCEDURE salsa_init_arrays
802    END INTERFACE salsa_init_arrays
803!
804!-- Writing of binary output for restart runs  !!! renaming?!
[3006]805    INTERFACE salsa_wrd_local
806       MODULE PROCEDURE salsa_wrd_local
807    END INTERFACE salsa_wrd_local
[2505]808!
809!-- Reading of NAMELIST parameters
810    INTERFACE salsa_parin
811       MODULE PROCEDURE salsa_parin
812    END INTERFACE salsa_parin
813!
814!-- Reading of parameters for restart runs
[3006]815    INTERFACE salsa_rrd_local
816       MODULE PROCEDURE salsa_rrd_local
817    END INTERFACE salsa_rrd_local
[2505]818!
819!-- Swapping of time levels (required for prognostic variables)
820    INTERFACE salsa_swap_timelevel
821       MODULE PROCEDURE salsa_swap_timelevel
822    END INTERFACE salsa_swap_timelevel
[3864]823!
824!-- Interface between PALM and salsa
[2505]825    INTERFACE salsa_driver
826       MODULE PROCEDURE salsa_driver
827    END INTERFACE salsa_driver
[3871]828
829!-- Actions salsa variables
830    INTERFACE salsa_actions
831       MODULE PROCEDURE salsa_actions
832       MODULE PROCEDURE salsa_actions_ij
833    END INTERFACE salsa_actions
[3864]834!
[3956]835!-- Non-advective processes (i.e. aerosol dynamic reactions) for salsa variables
836    INTERFACE salsa_non_advective_processes
837       MODULE PROCEDURE salsa_non_advective_processes
838       MODULE PROCEDURE salsa_non_advective_processes_ij
839    END INTERFACE salsa_non_advective_processes
840!
841!-- Exchange horiz for salsa variables
842    INTERFACE salsa_exchange_horiz_bounds
843       MODULE PROCEDURE salsa_exchange_horiz_bounds
844    END INTERFACE salsa_exchange_horiz_bounds
845!
[3864]846!-- Prognostics equations for salsa variables
847    INTERFACE salsa_prognostic_equations
848       MODULE PROCEDURE salsa_prognostic_equations
849       MODULE PROCEDURE salsa_prognostic_equations_ij
850    END INTERFACE salsa_prognostic_equations
851!
852!-- Tendency salsa variables
[2738]853    INTERFACE salsa_tendency
854       MODULE PROCEDURE salsa_tendency
[3064]855       MODULE PROCEDURE salsa_tendency_ij
[2738]856    END INTERFACE salsa_tendency
[3864]857
858
[2505]859    SAVE
860
861    PRIVATE
862!
[3864]863!-- Public functions:
864    PUBLIC salsa_boundary_conds, salsa_check_data_output, salsa_check_parameters,                  &
865           salsa_3d_data_averaging, salsa_data_output_2d, salsa_data_output_3d,                    &
866           salsa_data_output_mask, salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver,      &
867           salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin,        &
[3871]868           salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local,     &
[3956]869           salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds
[2505]870!
871!-- Public parameters, constants and initial values
[3864]872    PUBLIC bc_am_t_val, bc_an_t_val, bc_gt_t_val, dots_salsa, dt_salsa,                            &
873           ibc_salsa_b, last_salsa_time, lsdepo, nest_salsa, salsa, salsa_gases_from_chem,         &
874           skip_time_do_salsa
[2505]875!
876!-- Public prognostic variables
[3864]877    PUBLIC aerosol_mass, aerosol_number, gconc_2, mconc_2, nbins_aerosol, ncc, ncomponents_mass,   &
878           nclim, nconc_2, ngases_salsa, prtcl, ra_dry, salsa_gas, sedim_vd
[2505]879
[3864]880
[2505]881 CONTAINS
882
883!------------------------------------------------------------------------------!
884! Description:
885! ------------
886!> Parin for &salsa_par for new modules
887!------------------------------------------------------------------------------!
888 SUBROUTINE salsa_parin
889
890    IMPLICIT NONE
891
[3864]892    CHARACTER(LEN=80) ::  line   !< dummy string that contains the current line
893                                  !< of the parameter file
894
895    NAMELIST /salsa_parameters/      aerosol_flux_dpg, aerosol_flux_mass_fracs_a,                  &
896                                     aerosol_flux_mass_fracs_b, aerosol_flux_sigmag,               &
[4117]897                                     advect_particle_water, bc_salsa_b, bc_salsa_t,                &
898                                     decycle_salsa_lr, decycle_method_salsa, decycle_salsa_ns,     &
899                                     depo_pcm_par, depo_pcm_type,                                  &
[3864]900                                     depo_surf_par, dpg, dt_salsa, feedback_to_palm, h2so4_init,   &
[3899]901                                     hno3_init, init_gases_type, init_aerosol_type, listspec,      &
902                                     mass_fracs_a, mass_fracs_b, n_lognorm, nbin, nest_salsa, nf2a,&
903                                     nh3_init, nj3, nlcnd, nlcndgas, nlcndh2oae, nlcoag, nldepo,   &
904                                     nldepo_pcm,  nldepo_surf, nldistupdate, nsnucl, ocnv_init,    &
905                                     ocsv_init, read_restart_data_salsa, reglim, salsa,            &
906                                     salsa_emission_mode, sigmag, skip_time_do_salsa,              &
907                                     surface_aerosol_flux, van_der_waals_coagc, write_binary_salsa
[3864]908
[2505]909    line = ' '
910!
911!-- Try to find salsa package
912    REWIND ( 11 )
913    line = ' '
[3308]914    DO WHILE ( INDEX( line, '&salsa_parameters' ) == 0 )
[2505]915       READ ( 11, '(A)', END=10 )  line
916    ENDDO
917    BACKSPACE ( 11 )
918!
919!-- Read user-defined namelist
[3308]920    READ ( 11, salsa_parameters )
[2505]921!
[3582]922!-- Enable salsa (salsa switch in modules.f90)
[2505]923    salsa = .TRUE.
924
925 10 CONTINUE
[3864]926
[2505]927 END SUBROUTINE salsa_parin
928
929!------------------------------------------------------------------------------!
930! Description:
931! ------------
932!> Check parameters routine for salsa.
933!------------------------------------------------------------------------------!
934 SUBROUTINE salsa_check_parameters
935
[3864]936    USE control_parameters,                                                                        &
[4117]937        ONLY:  humidity
[3864]938
[2505]939    IMPLICIT NONE
[3864]940
[2505]941!
942!-- Checks go here (cf. check_parameters.f90).
943    IF ( salsa  .AND.  .NOT.  humidity )  THEN
[3864]944       WRITE( message_string, * ) 'salsa = ', salsa, ' is not allowed with humidity = ', humidity
945       CALL message( 'salsa_check_parameters', 'PA0594', 1, 2, 0, 6, 0 )
[2505]946    ENDIF
[3864]947
[2777]948    IF ( bc_salsa_b == 'dirichlet' )  THEN
[2505]949       ibc_salsa_b = 0
[2777]950    ELSEIF ( bc_salsa_b == 'neumann' )  THEN
[2505]951       ibc_salsa_b = 1
952    ELSE
[3864]953       message_string = 'unknown boundary condition: bc_salsa_b = "' // TRIM( bc_salsa_t ) // '"'
954       CALL message( 'salsa_check_parameters', 'PA0595', 1, 2, 0, 6, 0 )
[2505]955    ENDIF
[3864]956
[2505]957    IF ( bc_salsa_t == 'dirichlet' )  THEN
958       ibc_salsa_t = 0
959    ELSEIF ( bc_salsa_t == 'neumann' )  THEN
960       ibc_salsa_t = 1
[3864]961    ELSEIF ( bc_salsa_t == 'nested' )  THEN
962       ibc_salsa_t = 2
[2505]963    ELSE
[3864]964       message_string = 'unknown boundary condition: bc_salsa_t = "' // TRIM( bc_salsa_t ) // '"'
965       CALL message( 'salsa_check_parameters', 'PA0596', 1, 2, 0, 6, 0 )
[2505]966    ENDIF
[3864]967
[3331]968    IF ( nj3 < 1  .OR.  nj3 > 3 )  THEN
969       message_string = 'unknown nj3 (must be 1-3)'
[3864]970       CALL message( 'salsa_check_parameters', 'PA0597', 1, 2, 0, 6, 0 )
[3331]971    ENDIF
[3864]972
[3899]973    IF ( salsa_emission_mode /= 'no_emission'  .AND.  ibc_salsa_b  == 0 ) THEN
974       message_string = 'salsa_emission_mode /= "no_emission" requires bc_salsa_b = "Neumann"'
[3864]975       CALL message( 'salsa_check_parameters','PA0598', 1, 2, 0, 6, 0 )
976    ENDIF
977
[3899]978    IF ( salsa_emission_mode /= 'no_emission' )  include_emission = .TRUE.
979
[2505]980 END SUBROUTINE salsa_check_parameters
981
982!------------------------------------------------------------------------------!
983!
984! Description:
985! ------------
986!> Subroutine defining appropriate grid for netcdf variables.
[3864]987!> It is called out from subroutine netcdf.
[2505]988!> Same grid as for other scalars (see netcdf_interface_mod.f90)
989!------------------------------------------------------------------------------!
990 SUBROUTINE salsa_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
[3864]991
[2505]992    IMPLICIT NONE
993
[3864]994    CHARACTER(LEN=*), INTENT(OUT) ::  grid_x   !<
995    CHARACTER(LEN=*), INTENT(OUT) ::  grid_y   !<
996    CHARACTER(LEN=*), INTENT(OUT) ::  grid_z   !<
997    CHARACTER(LEN=*), INTENT(IN)  ::  var      !<
998
[2505]999    LOGICAL, INTENT(OUT) ::  found   !<
[3864]1000
[2505]1001    found  = .TRUE.
1002!
1003!-- Check for the grid
1004
[3412]1005    IF ( var(1:2) == 'g_' )  THEN
[3864]1006       grid_x = 'x'
1007       grid_y = 'y'
1008       grid_z = 'zu'
[3412]1009    ELSEIF ( var(1:4) == 'LDSA' )  THEN
[3864]1010       grid_x = 'x'
1011       grid_y = 'y'
[3412]1012       grid_z = 'zu'
1013    ELSEIF ( var(1:5) == 'm_bin' )  THEN
[3864]1014       grid_x = 'x'
1015       grid_y = 'y'
[3412]1016       grid_z = 'zu'
1017    ELSEIF ( var(1:5) == 'N_bin' )  THEN
[3864]1018       grid_x = 'x'
1019       grid_y = 'y'
[3412]1020       grid_z = 'zu'
1021    ELSEIF ( var(1:4) == 'Ntot' ) THEN
[3864]1022       grid_x = 'x'
1023       grid_y = 'y'
[3412]1024       grid_z = 'zu'
1025    ELSEIF ( var(1:2) == 'PM' )  THEN
[3864]1026       grid_x = 'x'
1027       grid_y = 'y'
[3412]1028       grid_z = 'zu'
1029    ELSEIF ( var(1:2) == 's_' )  THEN
[3864]1030       grid_x = 'x'
1031       grid_y = 'y'
[3412]1032       grid_z = 'zu'
1033    ELSE
1034       found  = .FALSE.
1035       grid_x = 'none'
1036       grid_y = 'none'
1037       grid_z = 'none'
1038    ENDIF
[2505]1039
1040 END SUBROUTINE salsa_define_netcdf_grid
1041
1042!------------------------------------------------------------------------------!
1043! Description:
1044! ------------
1045!> Header output for new module
1046!------------------------------------------------------------------------------!
1047 SUBROUTINE salsa_header( io )
1048
[4012]1049    USE indices,                                                                                   &
1050        ONLY:  nx, ny, nz
1051
[2505]1052    IMPLICIT NONE
1053 
1054    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
1055!
1056!-- Write SALSA header
1057    WRITE( io, 1 )
1058    WRITE( io, 2 ) skip_time_do_salsa
1059    WRITE( io, 3 ) dt_salsa
[4012]1060    WRITE( io, 4 )  nz, ny, nx, nbins_aerosol
[2869]1061    IF ( advect_particle_water )  THEN
[4012]1062       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol,              &
[2869]1063                        advect_particle_water
1064    ELSE
[3864]1065       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncc*nbins_aerosol, advect_particle_water
[2869]1066    ENDIF
1067    IF ( .NOT. salsa_gases_from_chem )  THEN
[3864]1068       WRITE( io, 6 )  SHAPE( aerosol_mass(1)%conc ), ngases_salsa, salsa_gases_from_chem
[2869]1069    ENDIF
[3864]1070    WRITE( io, 7 )
[3899]1071    IF ( nsnucl > 0 )   WRITE( io, 8 ) nsnucl, nj3
1072    IF ( nlcoag )       WRITE( io, 9 )
1073    IF ( nlcnd )        WRITE( io, 10 ) nlcndgas, nlcndh2oae
1074    IF ( lspartition )  WRITE( io, 11 )
1075    IF ( nldepo )       WRITE( io, 12 ) nldepo_pcm, nldepo_surf
[3864]1076    WRITE( io, 13 )  reglim, nbin, bin_low_limits
[3899]1077    IF ( init_aerosol_type == 0 )  WRITE( io, 14 ) nsect
[3864]1078    WRITE( io, 15 ) ncc, listspec, mass_fracs_a, mass_fracs_b
[3005]1079    IF ( .NOT. salsa_gases_from_chem )  THEN
[3864]1080       WRITE( io, 16 ) ngases_salsa, h2so4_init, hno3_init, nh3_init, ocnv_init, ocsv_init
[3005]1081    ENDIF
[3899]1082    WRITE( io, 17 )  init_aerosol_type, init_gases_type
1083    IF ( init_aerosol_type == 0 )  THEN
[3864]1084       WRITE( io, 18 )  dpg, sigmag, n_lognorm
[2505]1085    ELSE
[3864]1086       WRITE( io, 19 )
[2505]1087    ENDIF
[3864]1088    IF ( nest_salsa )  WRITE( io, 20 )  nest_salsa
1089    WRITE( io, 21 ) salsa_emission_mode
[3899]1090    IF ( salsa_emission_mode == 'uniform' )  THEN
1091       WRITE( io, 22 ) surface_aerosol_flux, aerosol_flux_dpg, aerosol_flux_sigmag,                &
1092                       aerosol_flux_mass_fracs_a
1093    ENDIF
1094    IF ( SUM( aerosol_flux_mass_fracs_b ) > 0.0_wp  .OR. salsa_emission_mode == 'read_from_file' ) &
1095    THEN
1096       WRITE( io, 23 )
1097    ENDIF
[2505]1098
[3864]10991   FORMAT (//' SALSA information:'/                                                               &
[2505]1100              ' ------------------------------'/)
[2561]11012   FORMAT   ('    Starts at: skip_time_do_salsa = ', F10.2, '  s')
[2505]11023   FORMAT  (/'    Timestep: dt_salsa = ', F6.2, '  s')
[3864]11034   FORMAT  (/'    Array shape (z,y,x,bins):'/                                                     &
[2869]1104              '       aerosol_number:  ', 4(I3)) 
[3864]11055   FORMAT  (/'       aerosol_mass:    ', 4(I3),/                                                  &
[3005]1106              '       (advect_particle_water = ', L1, ')')
[3864]11076   FORMAT   ('       salsa_gas: ', 4(I3),/                                                        &
[3005]1108              '       (salsa_gases_from_chem = ', L1, ')')
[3864]11097   FORMAT  (/'    Aerosol dynamic processes included: ')
11108   FORMAT  (/'       nucleation (scheme = ', I1, ' and J3 parametrization = ', I1, ')')
11119   FORMAT  (/'       coagulation')
111210  FORMAT  (/'       condensation (of precursor gases = ', L1, ' and water vapour = ', L1, ')' )
111311  FORMAT  (/'       dissolutional growth by HNO3 and NH3')
111412  FORMAT  (/'       dry deposition (on vegetation = ', L1, ' and on topography = ', L1, ')')
111513  FORMAT  (/'    Aerosol bin subrange limits (in metres): ',  3(ES10.2E3), /                     &
1116              '    Number of size bins for each aerosol subrange: ', 2I3,/                         &
[3780]1117              '    Aerosol bin limits (in metres): ', 9(ES10.2E3))
[3864]111814  FORMAT   ('    Initial number concentration in bins at the lowest level (#/m**3):', 9(ES10.2E3))
111915  FORMAT  (/'    Number of chemical components used: ', I1,/                                     &
1120              '       Species: ',7(A6),/                                                           &
1121              '    Initial relative contribution of each species to particle volume in:',/         &
1122              '       a-bins: ', 7(F6.3),/                                                         &
[3005]1123              '       b-bins: ', 7(F6.3))
[3864]112416  FORMAT  (/'    Number of gaseous tracers used: ', I1,/                                         &
1125              '    Initial gas concentrations:',/                                                  &
1126              '       H2SO4: ',ES12.4E3, ' #/m**3',/                                               &
1127              '       HNO3:  ',ES12.4E3, ' #/m**3',/                                               &
1128              '       NH3:   ',ES12.4E3, ' #/m**3',/                                               &
1129              '       OCNV:  ',ES12.4E3, ' #/m**3',/                                               &
[2505]1130              '       OCSV:  ',ES12.4E3, ' #/m**3')
[3864]113117   FORMAT (/'   Initialising concentrations: ', /                                                &
[3899]1132              '      Aerosol size distribution: init_aerosol_type = ', I1,/                        &
1133              '      Gas concentrations: init_gases_type = ', I1 )
[3864]113418   FORMAT ( '      Mode diametres: dpg(nmod) = ', 7(F7.3), ' (m)', /                             &
1135              '      Standard deviation: sigmag(nmod) = ', 7(F7.2),/                               &
1136              '      Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3), ' (#/m3)' )
113719   FORMAT (/'      Size distribution read from a file.')
113820   FORMAT (/'   Nesting for salsa variables: ', L1 )
113921   FORMAT (/'   Emissions: salsa_emission_mode = ', A )
[3899]114022   FORMAT (/'      surface_aerosol_flux = ', ES12.4E3, ' #/m**2/s', /                            &
1141              '      aerosol_flux_dpg     =  ', 7(F7.3), ' (m)', /                                 &
1142              '      aerosol_flux_sigmag  =  ', 7(F7.2), /                                         &
1143              '      aerosol_mass_fracs_a =  ', 7(ES12.4E3) )
114423   FORMAT (/'      (currently all emissions are soluble!)')
[2505]1145
1146 END SUBROUTINE salsa_header
1147
1148!------------------------------------------------------------------------------!
1149! Description:
1150! ------------
1151!> Allocate SALSA arrays and define pointers if required
1152!------------------------------------------------------------------------------!
1153 SUBROUTINE salsa_init_arrays
1154
[4117]1155    USE advec_ws,                                                                                  &
1156        ONLY: ws_init_flags_scalar
[3864]1157
1158    USE surface_mod,                                                                               &
1159        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
1160
[2505]1161    IMPLICIT NONE
[3864]1162
1163    INTEGER(iwp) ::  gases_available !< Number of available gas components in the chemistry model
1164    INTEGER(iwp) ::  i               !< loop index for allocating
1165    INTEGER(iwp) ::  l               !< loop index for allocating: surfaces
1166    INTEGER(iwp) ::  lsp             !< loop index for chem species in the chemistry model
1167
[2754]1168    gases_available = 0
[2505]1169!
1170!-- Allocate prognostic variables (see salsa_swap_timelevel)
[3864]1171!
1172!-- Set derived indices:
1173!-- (This does the same as the subroutine salsa_initialize in SALSA/UCLALES-SALSA)
1174    start_subrange_1a = 1  ! 1st index of subrange 1a
1175    start_subrange_2a = start_subrange_1a + nbin(1)  ! 1st index of subrange 2a
1176    end_subrange_1a   = start_subrange_2a - 1        ! last index of subrange 1a
1177    end_subrange_2a   = end_subrange_1a + nbin(2)    ! last index of subrange 2a
[3636]1178
[2505]1179!
[3864]1180!-- If the fraction of insoluble aerosols in subrange 2 is zero: do not allocate arrays for them
[2777]1181    IF ( nf2a > 0.999999_wp  .AND.  SUM( mass_fracs_b ) < 0.00001_wp )  THEN
[2505]1182       no_insoluble = .TRUE.
[3864]1183       start_subrange_2b = end_subrange_2a+1  ! 1st index of subrange 2b
1184       end_subrange_2b   = end_subrange_2a    ! last index of subrange 2b
[2505]1185    ELSE
[3864]1186       start_subrange_2b = start_subrange_2a + nbin(2)  ! 1st index of subrange 2b
1187       end_subrange_2b   = end_subrange_2a + nbin(2)    ! last index of subrange 2b
[2505]1188    ENDIF
[3864]1189
1190    nbins_aerosol = end_subrange_2b   ! total number of aerosol size bins
1191!
[2505]1192!-- Create index tables for different aerosol components
[2754]1193    CALL component_index_constructor( prtcl, ncc, maxspec, listspec )
[3864]1194
1195    ncomponents_mass = ncc
1196    IF ( advect_particle_water )  ncomponents_mass = ncc + 1  ! Add water
1197
[2505]1198!
1199!-- Allocate:
[3864]1200    ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass),                    &
[3899]1201              bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),&
[3864]1202              nsect(nbins_aerosol), massacc(nbins_aerosol) )
1203    ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) )
1204    IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1205    ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
1206
1207!
[2505]1208!-- Aerosol number concentration
[3864]1209    ALLOCATE( aerosol_number(nbins_aerosol) )
1210    ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
1211              nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
1212              nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
[2525]1213    nconc_1 = 0.0_wp
1214    nconc_2 = 0.0_wp
1215    nconc_3 = 0.0_wp
[3864]1216
1217    DO i = 1, nbins_aerosol
[2505]1218       aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => nconc_1(:,:,:,i)
1219       aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => nconc_2(:,:,:,i)
1220       aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i)
[3899]1221       ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                         &
1222                 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                         &
1223                 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
1224                 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
1225                 aerosol_number(i)%init(nzb:nzt+1),                                                &
[2505]1226                 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
[4012]1227       aerosol_number(i)%init = nclim
[3899]1228       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
1229          ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
1230       ENDIF
[3864]1231    ENDDO
1232
1233!
1234!-- Aerosol mass concentration
1235    ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) )
1236    ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
1237              mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
1238              mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) )
[2525]1239    mconc_1 = 0.0_wp
1240    mconc_2 = 0.0_wp
1241    mconc_3 = 0.0_wp
[3864]1242
1243    DO i = 1, ncomponents_mass*nbins_aerosol
[2505]1244       aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => mconc_1(:,:,:,i)
1245       aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => mconc_2(:,:,:,i)
[3864]1246       aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i)
1247       ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                           &
1248                 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                           &
1249                 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
1250                 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
1251                 aerosol_mass(i)%init(nzb:nzt+1),                                                  &
[2505]1252                 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
[4012]1253       aerosol_mass(i)%init = mclim
[3899]1254       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
1255          ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
1256       ENDIF
[2505]1257    ENDDO
[3864]1258
[2505]1259!
[3864]1260!-- Surface fluxes: answs = aerosol number, amsws = aerosol mass
[2505]1261!
1262!-- Horizontal surfaces: default type
1263    DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
[3864]1264       ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) )
1265       ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) )
[2505]1266       surf_def_h(l)%answs = 0.0_wp
1267       surf_def_h(l)%amsws = 0.0_wp
1268    ENDDO
1269!
[3864]1270!-- Horizontal surfaces: natural type
1271    ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) )
1272    ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) )
1273    surf_lsm_h%answs = 0.0_wp
1274    surf_lsm_h%amsws = 0.0_wp
1275!
1276!-- Horizontal surfaces: urban type
1277    ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) )
1278    ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) )
1279    surf_usm_h%answs = 0.0_wp
1280    surf_usm_h%amsws = 0.0_wp
1281
1282!
1283!-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing
1284    DO  l = 0, 3
1285       ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) )
[2505]1286       surf_def_v(l)%answs = 0.0_wp
[3864]1287       ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
[2505]1288       surf_def_v(l)%amsws = 0.0_wp
[3864]1289
1290       ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) )
1291       surf_lsm_v(l)%answs = 0.0_wp
1292       ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1293       surf_lsm_v(l)%amsws = 0.0_wp
1294
1295       ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) )
1296       surf_usm_v(l)%answs = 0.0_wp
1297       ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
1298       surf_usm_v(l)%amsws = 0.0_wp
1299
1300    ENDDO
1301
[2754]1302!
1303!-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV)
1304!-- (number concentration (#/m3) )
1305!
1306!-- If chemistry is on, read gas phase concentrations from there. Otherwise,
[3308]1307!-- allocate salsa_gas array.
[2505]1308
[3864]1309    IF ( air_chemistry )  THEN
[2754]1310       DO  lsp = 1, nvar
[3864]1311          SELECT CASE ( TRIM( chem_species(lsp)%name ) )
1312             CASE ( 'H2SO4', 'h2so4' )
1313                gases_available = gases_available + 1
1314                gas_index_chem(1) = lsp
1315             CASE ( 'HNO3', 'hno3' )
1316                gases_available = gases_available + 1
1317                gas_index_chem(2) = lsp
1318             CASE ( 'NH3', 'nh3' )
1319                gases_available = gases_available + 1
1320                gas_index_chem(3) = lsp
1321             CASE ( 'OCNV', 'ocnv' )
1322                gases_available = gases_available + 1
1323                gas_index_chem(4) = lsp
1324             CASE ( 'OCSV', 'ocsv' )
1325                gases_available = gases_available + 1
1326                gas_index_chem(5) = lsp
1327          END SELECT
[2754]1328       ENDDO
1329
[3864]1330       IF ( gases_available == ngases_salsa )  THEN
[2754]1331          salsa_gases_from_chem = .TRUE.
1332       ELSE
[3864]1333          WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// &
1334                                     'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)'
1335       CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 )
[2754]1336       ENDIF
[2505]1337
[2754]1338    ELSE
1339
[3864]1340       ALLOCATE( salsa_gas(ngases_salsa) )
1341       ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
1342                 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
1343                 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
[2754]1344       gconc_1 = 0.0_wp
1345       gconc_2 = 0.0_wp
1346       gconc_3 = 0.0_wp
[3864]1347
1348       DO i = 1, ngases_salsa
[3308]1349          salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => gconc_1(:,:,:,i)
1350          salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => gconc_2(:,:,:,i)
1351          salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i)
1352          ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
1353                    salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
1354                    salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1355                    salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1356                    salsa_gas(i)%init(nzb:nzt+1),                              &
1357                    salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
[4012]1358          salsa_gas(i)%init = nclim
[3899]1359          IF ( include_emission )  ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
[3864]1360       ENDDO
[2754]1361!
1362!--    Surface fluxes: gtsws = gaseous tracer flux
1363!
1364!--    Horizontal surfaces: default type
1365       DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
[3864]1366          ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) )
[2754]1367          surf_def_h(l)%gtsws = 0.0_wp
1368       ENDDO
[3864]1369!--    Horizontal surfaces: natural type
1370       ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) )
1371       surf_lsm_h%gtsws = 0.0_wp
1372!--    Horizontal surfaces: urban type
1373       ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) )
1374       surf_usm_h%gtsws = 0.0_wp
[2754]1375!
1376!--    Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
1377!--    westward (l=3) facing
[3864]1378       DO  l = 0, 3
1379          ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) )
[2754]1380          surf_def_v(l)%gtsws = 0.0_wp
[3864]1381          ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) )
1382          surf_lsm_v(l)%gtsws = 0.0_wp
1383          ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) )
1384          surf_usm_v(l)%gtsws = 0.0_wp
[2754]1385       ENDDO
1386    ENDIF
[3864]1387
[3871]1388    IF ( ws_scheme_sca )  THEN
1389
1390       IF ( salsa )  THEN
1391          ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1392          sums_salsa_ws_l = 0.0_wp
1393       ENDIF
1394
1395    ENDIF
[4117]1396!
1397!-- Set control flags for decycling only at lateral boundary cores. Within the inner cores the
1398!-- decycle flag is set to .FALSE.. Even though it does not affect the setting of chemistry boundary
1399!-- conditions, this flag is used to set advection control flags appropriately.
1400    decycle_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0  .OR.  nxr == nx )
1401    decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0  .OR.  nyn == ny )
1402!
1403!-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have
1404!-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce
1405!-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
1406!-- applied near these boundaries. To get rid-off this, set-up additional flags that control the
1407!-- order of the scalar advection scheme near the lateral boundaries for passive scalars with
1408!-- decycling.
1409    IF ( scalar_advec == 'ws-scheme' )  THEN
1410       ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1411!
1412!--    In case of decycling, set Neuman boundary conditions for wall_flags_0 bit 31 instead of
1413!--    cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see
1414!--    the following comment). Note, since several also other modules may access this bit but may
1415!--    have other boundary conditions, the original value of wall_flags_0 bit 31 must not be
1416!--    modified. Hence, store the boundary conditions directly on salsa_advc_flags_s.
1417!--    salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used
1418!--    to control the numerical order.
1419!--    Initialize with flag 31 only.
1420       salsa_advc_flags_s = 0
1421       salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_0, 31 ) )
[3871]1422
[4117]1423       IF ( decycle_salsa_ns )  THEN
1424          IF ( nys == 0 )  THEN
1425             DO  i = 1, nbgp
1426                salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ),   &
1427                                                       IBCLR( salsa_advc_flags_s(:,nys,:), 31 ),   &
1428                                                       BTEST( salsa_advc_flags_s(:,nys,:), 31 ) )
1429             ENDDO
1430          ENDIF
1431          IF ( nyn == ny )  THEN
1432             DO  i = 1, nbgp
1433                salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ),   &
1434                                                       IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ),   &
1435                                                       BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) )
1436             ENDDO
1437          ENDIF
1438       ENDIF
1439       IF ( decycle_salsa_lr )  THEN
1440          IF ( nxl == 0 )  THEN
1441             DO  i = 1, nbgp
1442                salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ),   &
1443                                                       IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ),   &
1444                                                       BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) )
1445             ENDDO
1446          ENDIF
1447          IF ( nxr == nx )  THEN
1448             DO  i = 1, nbgp
1449                salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ),   &
1450                                                       IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ),   &
1451                                                       BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) )
1452             ENDDO
1453          ENDIF
1454       ENDIF
1455!
1456!--    To initialise the advection flags appropriately, pass the boundary flags to
1457!--    ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive
1458!--    scalar is being treated and the horizontal advection terms are degraded already 2 grid points
1459!--    before the lateral boundary. Also, extended degradation zones are applied, where
1460!--    horizontal advection of scalars is discretised by the first-order scheme at all grid points
1461!--    in the vicinity of buildings (<= 3 grid points). Even though no building is within the
1462!--    numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order
1463!--    of the horizontal advection scheme is successively upgraded.
1464!--    These degradations of the advection scheme are done to avoid stationary numerical
1465!--    oscillations, which are responsible for high concentration maxima that may appear e.g. under
1466!--    shear-free stable conditions.
1467       CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_salsa_lr,    &
1468                                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_salsa_ns,    &
1469                                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_salsa_lr,    &
1470                                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_salsa_ns,    &
1471                                  salsa_advc_flags_s, .TRUE. )
1472    ENDIF
1473
1474
[2505]1475 END SUBROUTINE salsa_init_arrays
1476
1477!------------------------------------------------------------------------------!
1478! Description:
1479! ------------
[3864]1480!> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA.
[2505]1481!> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are
1482!> also merged here.
1483!------------------------------------------------------------------------------!
1484 SUBROUTINE salsa_init
1485
1486    IMPLICIT NONE
[3864]1487
1488    INTEGER(iwp) :: i   !<
1489    INTEGER(iwp) :: ib  !< loop index for aerosol number bins
1490    INTEGER(iwp) :: ic  !< loop index for aerosol mass bins
1491    INTEGER(iwp) :: ig  !< loop index for gases
1492    INTEGER(iwp) :: ii  !< index for indexing
1493    INTEGER(iwp) :: j   !<
1494
[3885]1495    IF ( debug_output )  CALL debug_message( 'salsa_init', 'start' )
[3864]1496
[2734]1497    bin_low_limits = 0.0_wp
[3864]1498    k_topo_top     = 0
[3125]1499    nsect          = 0.0_wp
[3864]1500    massacc        = 1.0_wp
1501
[2505]1502!
[2525]1503!-- Indices for chemical components used (-1 = not used)
[3864]1504    ii = 0
[2525]1505    IF ( is_used( prtcl, 'SO4' ) )  THEN
[3864]1506       index_so4 = get_index( prtcl,'SO4' )
1507       ii = ii + 1
[2525]1508    ENDIF
1509    IF ( is_used( prtcl,'OC' ) )  THEN
[3864]1510       index_oc = get_index(prtcl, 'OC')
1511       ii = ii + 1
[2525]1512    ENDIF
1513    IF ( is_used( prtcl, 'BC' ) )  THEN
[3864]1514       index_bc = get_index( prtcl, 'BC' )
1515       ii = ii + 1
[2525]1516    ENDIF
1517    IF ( is_used( prtcl, 'DU' ) )  THEN
[3864]1518       index_du = get_index( prtcl, 'DU' )
1519       ii = ii + 1
[2525]1520    ENDIF
1521    IF ( is_used( prtcl, 'SS' ) )  THEN
[3864]1522       index_ss = get_index( prtcl, 'SS' )
1523       ii = ii + 1
[2525]1524    ENDIF
1525    IF ( is_used( prtcl, 'NO' ) )  THEN
[3864]1526       index_no = get_index( prtcl, 'NO' )
1527       ii = ii + 1
[2525]1528    ENDIF
1529    IF ( is_used( prtcl, 'NH' ) )  THEN
[3864]1530       index_nh = get_index( prtcl, 'NH' )
1531       ii = ii + 1
[2525]1532    ENDIF
[3864]1533!
[2525]1534!-- All species must be known
[3864]1535    IF ( ii /= ncc )  THEN
1536       message_string = 'Unknown aerosol species/component(s) given in the initialization'
1537       CALL message( 'salsa_mod: salsa_init', 'PA0600', 1, 2, 0, 6, 0 )
[2525]1538    ENDIF
1539!
[3864]1540!-- Partition and dissolutional growth by gaseous HNO3 and NH3
1541    IF ( index_no > 0  .AND.  index_nh > 0  .AND.  index_so4 > 0 )  lspartition = .TRUE.
1542!
[2505]1543!-- Initialise
[2525]1544!
1545!-- Aerosol size distribution (TYPE t_section)
[2744]1546    aero(:)%dwet     = 1.0E-10_wp
1547    aero(:)%veqh2o   = 1.0E-10_wp
1548    aero(:)%numc     = nclim
1549    aero(:)%core     = 1.0E-10_wp
[3864]1550    DO ic = 1, maxspec+1    ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
1551       aero(:)%volc(ic) = 0.0_wp
[2505]1552    ENDDO
[3864]1553
1554    IF ( nldepo )  sedim_vd = 0.0_wp
1555
[3582]1556    IF ( .NOT. salsa_gases_from_chem )  THEN
[4012]1557       IF ( .NOT. read_restart_data_salsa )  THEN
1558          salsa_gas(1)%conc = h2so4_init
1559          salsa_gas(2)%conc = hno3_init
1560          salsa_gas(3)%conc = nh3_init
1561          salsa_gas(4)%conc = ocnv_init
1562          salsa_gas(5)%conc = ocsv_init
1563       ENDIF
[3864]1564       DO  ig = 1, ngases_salsa
1565          salsa_gas(ig)%conc_p    = 0.0_wp
1566          salsa_gas(ig)%tconc_m   = 0.0_wp
1567          salsa_gas(ig)%flux_s    = 0.0_wp
1568          salsa_gas(ig)%diss_s    = 0.0_wp
1569          salsa_gas(ig)%flux_l    = 0.0_wp
1570          salsa_gas(ig)%diss_l    = 0.0_wp
1571          salsa_gas(ig)%sums_ws_l = 0.0_wp
[4012]1572          salsa_gas(ig)%conc_p    = salsa_gas(ig)%conc
[2525]1573       ENDDO
[2505]1574!
[4012]1575!--    Set initial value for gas compound tracer
[3864]1576       salsa_gas(1)%init = h2so4_init
1577       salsa_gas(2)%init = hno3_init
1578       salsa_gas(3)%init = nh3_init
1579       salsa_gas(4)%init = ocnv_init
1580       salsa_gas(5)%init = ocsv_init
[3582]1581    ENDIF
1582!
1583!-- Aerosol radius in each bin: dry and wet (m)
[3864]1584    ra_dry = 1.0E-10_wp
1585!
1586!-- Initialise aerosol tracers
[3582]1587    aero(:)%vhilim   = 0.0_wp
1588    aero(:)%vlolim   = 0.0_wp
1589    aero(:)%vratiohi = 0.0_wp
1590    aero(:)%vratiolo = 0.0_wp
1591    aero(:)%dmid     = 0.0_wp
[2607]1592!
[3582]1593!-- Initialise the sectional particle size distribution
[3864]1594    CALL set_sizebins
[2505]1595!
[3864]1596!-- Initialise location-dependent aerosol size distributions and chemical compositions:
1597    CALL aerosol_init
[4012]1598
[3864]1599!-- Initalisation run of SALSA + calculate the vertical top index of the topography
[3582]1600    DO  i = nxl, nxr
1601       DO  j = nys, nyn
[3864]1602
1603          k_topo_top(j,i) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,j,i), 12 ) ), DIM = 1 ) - 1
1604
[3582]1605          CALL salsa_driver( i, j, 1 )
1606          CALL salsa_diagnostics( i, j )
1607       ENDDO
[3864]1608    ENDDO
[4012]1609
1610    DO  ib = 1, nbins_aerosol
1611       aerosol_number(ib)%conc_p    = aerosol_number(ib)%conc
1612       aerosol_number(ib)%tconc_m   = 0.0_wp
1613       aerosol_number(ib)%flux_s    = 0.0_wp
1614       aerosol_number(ib)%diss_s    = 0.0_wp
1615       aerosol_number(ib)%flux_l    = 0.0_wp
1616       aerosol_number(ib)%diss_l    = 0.0_wp
1617       aerosol_number(ib)%sums_ws_l = 0.0_wp
1618    ENDDO
1619    DO  ic = 1, ncomponents_mass*nbins_aerosol
1620       aerosol_mass(ic)%conc_p    = aerosol_mass(ic)%conc
1621       aerosol_mass(ic)%tconc_m   = 0.0_wp
1622       aerosol_mass(ic)%flux_s    = 0.0_wp
1623       aerosol_mass(ic)%diss_s    = 0.0_wp
1624       aerosol_mass(ic)%flux_l    = 0.0_wp
1625       aerosol_mass(ic)%diss_l    = 0.0_wp
1626       aerosol_mass(ic)%sums_ws_l = 0.0_wp
1627    ENDDO
[2525]1628!
[4012]1629!
[3864]1630!-- Initialise the deposition scheme and surface types
1631    IF ( nldepo )  CALL init_deposition
1632
[3899]1633    IF ( include_emission )  THEN
[3864]1634!
1635!--    Read in and initialize emissions
1636       CALL salsa_emission_setup( .TRUE. )
[3924]1637       IF ( .NOT. salsa_gases_from_chem  .AND.  salsa_emission_mode == 'read_from_file' )  THEN
[3864]1638          CALL salsa_gas_emission_setup( .TRUE. )
1639       ENDIF
[2533]1640    ENDIF
[3864]1641
[3885]1642    IF ( debug_output )  CALL debug_message( 'salsa_init', 'end' )
[3864]1643
[2505]1644 END SUBROUTINE salsa_init
1645
1646!------------------------------------------------------------------------------!
1647! Description:
1648! ------------
1649!> Initializes particle size distribution grid by calculating size bin limits
1650!> and mid-size for *dry* particles in each bin. Called from salsa_initialize
1651!> (only at the beginning of simulation).
1652!> Size distribution described using:
[2744]1653!>   1) moving center method (subranges 1 and 2)
[2505]1654!>      (Jacobson, Atmos. Env., 31, 131-144, 1997)
[2744]1655!>   2) fixed sectional method (subrange 3)
1656!> Size bins in each subrange are spaced logarithmically
1657!> based on given subrange size limits and bin number.
[2505]1658!
1659!> Mona changed 06/2017: Use geometric mean diameter to describe the mean
1660!> particle diameter in a size bin, not the arithmeric mean which clearly
1661!> overestimates the total particle volume concentration.
1662!
1663!> Coded by:
1664!> Hannele Korhonen (FMI) 2005
1665!> Harri Kokkola (FMI) 2006
1666!
1667!> Bug fixes for box model + updated for the new aerosol datatype:
1668!> Juha Tonttila (FMI) 2014
1669!------------------------------------------------------------------------------!
1670 SUBROUTINE set_sizebins
[3864]1671
[2505]1672    IMPLICIT NONE
[3864]1673
1674    INTEGER(iwp) ::  cc  !< running index
1675    INTEGER(iwp) ::  dd  !< running index
1676
1677    REAL(wp) ::  ratio_d  !< ratio of the upper and lower diameter of subranges
[2505]1678!
[3864]1679!-- vlolim&vhilim: min & max *dry* volumes [fxm]
[2505]1680!-- dmid: bin mid *dry* diameter (m)
1681!-- vratiolo&vratiohi: volume ratio between the center and low/high limit
1682!
[2744]1683!-- 1) Size subrange 1:
[2505]1684    ratio_d = reglim(2) / reglim(1)   ! section spacing (m)
[3864]1685    DO  cc = start_subrange_1a, end_subrange_1a
1686       aero(cc)%vlolim = api6 * ( reglim(1) * ratio_d**( REAL( cc-1 ) / nbin(1) ) )**3
1687       aero(cc)%vhilim = api6 * ( reglim(1) * ratio_d**( REAL( cc ) / nbin(1) ) )**3
1688       aero(cc)%dmid = SQRT( ( aero(cc)%vhilim / api6 )**0.33333333_wp *                           &
1689                             ( aero(cc)%vlolim / api6 )**0.33333333_wp )
1690       aero(cc)%vratiohi = aero(cc)%vhilim / ( api6 * aero(cc)%dmid**3 )
1691       aero(cc)%vratiolo = aero(cc)%vlolim / ( api6 * aero(cc)%dmid**3 )
[2505]1692    ENDDO
1693!
[3864]1694!-- 2) Size subrange 2:
[2744]1695!-- 2.1) Sub-subrange 2a: high hygroscopicity
[2505]1696    ratio_d = reglim(3) / reglim(2)   ! section spacing
[3864]1697    DO  dd = start_subrange_2a, end_subrange_2a
1698       cc = dd - start_subrange_2a
1699       aero(dd)%vlolim = api6 * ( reglim(2) * ratio_d**( REAL( cc ) / nbin(2) ) )**3
1700       aero(dd)%vhilim = api6 * ( reglim(2) * ratio_d**( REAL( cc+1 ) / nbin(2) ) )**3
1701       aero(dd)%dmid = SQRT( ( aero(dd)%vhilim / api6 )**0.33333333_wp *                           &
1702                             ( aero(dd)%vlolim / api6 )**0.33333333_wp )
1703       aero(dd)%vratiohi = aero(dd)%vhilim / ( api6 * aero(dd)%dmid**3 )
1704       aero(dd)%vratiolo = aero(dd)%vlolim / ( api6 * aero(dd)%dmid**3 )
[2505]1705    ENDDO
[3864]1706!
[2744]1707!-- 2.2) Sub-subrange 2b: low hygroscopicity
[2505]1708    IF ( .NOT. no_insoluble )  THEN
[3864]1709       aero(start_subrange_2b:end_subrange_2b)%vlolim   = aero(start_subrange_2a:end_subrange_2a)%vlolim
1710       aero(start_subrange_2b:end_subrange_2b)%vhilim   = aero(start_subrange_2a:end_subrange_2a)%vhilim
1711       aero(start_subrange_2b:end_subrange_2b)%dmid     = aero(start_subrange_2a:end_subrange_2a)%dmid
1712       aero(start_subrange_2b:end_subrange_2b)%vratiohi = aero(start_subrange_2a:end_subrange_2a)%vratiohi
1713       aero(start_subrange_2b:end_subrange_2b)%vratiolo = aero(start_subrange_2a:end_subrange_2a)%vratiolo
[2505]1714    ENDIF
[3864]1715!
1716!-- Initialize the wet diameter with the bin dry diameter to avoid numerical problems later
[2505]1717    aero(:)%dwet = aero(:)%dmid
[2607]1718!
[3864]1719!-- Save bin limits (lower diameter) to be delivered to PALM if needed
1720    DO cc = 1, nbins_aerosol
1721       bin_low_limits(cc) = ( aero(cc)%vlolim / api6 )**0.33333333_wp
1722    ENDDO
1723
[2505]1724 END SUBROUTINE set_sizebins
[3864]1725
[2505]1726!------------------------------------------------------------------------------!
1727! Description:
1728! ------------
1729!> Initilize altitude-dependent aerosol size distributions and compositions.
1730!>
1731!> Mona added 06/2017: Correct the number and mass concentrations by normalizing
1732!< by the given total number and mass concentration.
1733!>
1734!> Tomi Raatikainen, FMI, 29.2.2016
1735!------------------------------------------------------------------------------!
1736 SUBROUTINE aerosol_init
[3864]1737
1738    USE netcdf_data_input_mod,                                                                     &
[4012]1739        ONLY:  close_input_file, get_attribute, get_variable,                                      &
1740               netcdf_data_input_get_dimension_length, open_read_file
[3864]1741
[2505]1742    IMPLICIT NONE
[2777]1743
[3864]1744    CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name  !< chemical component name
1745
1746    INTEGER(iwp) ::  ee        !< index: end
1747    INTEGER(iwp) ::  i         !< loop index: x-direction
1748    INTEGER(iwp) ::  ib        !< loop index: size bins
1749    INTEGER(iwp) ::  ic        !< loop index: chemical components
1750    INTEGER(iwp) ::  id_dyn    !< NetCDF id of PIDS_DYNAMIC_SALSA
1751    INTEGER(iwp) ::  ig        !< loop index: gases
1752    INTEGER(iwp) ::  j         !< loop index: y-direction
1753    INTEGER(iwp) ::  k         !< loop index: z-direction
1754    INTEGER(iwp) ::  lod_aero  !< level of detail of inital aerosol concentrations
1755    INTEGER(iwp) ::  pr_nbins  !< Number of aerosol size bins in file
1756    INTEGER(iwp) ::  pr_ncc    !< Number of aerosol chemical components in file
1757    INTEGER(iwp) ::  pr_nz     !< Number of vertical grid-points in file
1758    INTEGER(iwp) ::  prunmode  !< running mode of SALSA
1759    INTEGER(iwp) ::  ss        !< index: start
1760
1761    INTEGER(iwp), DIMENSION(maxspec) ::  cc_input_to_model
1762
1763    LOGICAL  ::  netcdf_extend = .FALSE. !< Flag: netcdf file exists
1764
1765    REAL(wp) ::  flag  !< flag to mask topography grid points
1766
1767    REAL(wp), DIMENSION(nbins_aerosol) ::  core   !< size of the bin mid aerosol particle
1768
1769    REAL(wp), DIMENSION(0:nz+1) ::  pnf2a   !< number fraction in 2a
1770    REAL(wp), DIMENSION(0:nz+1) ::  pmfoc1a !< mass fraction of OC in 1a
1771
1772    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol)   ::  pndist  !< size dist as a function of height (#/m3)
1773    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2a   !< mass distributions in subrange 2a
1774    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2b   !< mass distributions in subrange 2b
1775
1776    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_dmid  !< vertical profile of aerosol bin diameters
1777    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_z     !< z levels of profiles
1778
1779    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_a  !< mass fraction: a
1780    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_b  !< and b
1781
1782    cc_input_to_model = 0
[2744]1783    prunmode = 1
[2505]1784!
1785!-- Bin mean aerosol particle volume (m3)
1786    core(:) = 0.0_wp
[3864]1787    core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3
1788!
1789!-- Set concentrations to zero
[2505]1790    pndist(:,:)  = 0.0_wp
[3864]1791    pnf2a(:)     = nf2a
1792    pmf2a(:,:)   = 0.0_wp
1793    pmf2b(:,:)   = 0.0_wp
1794    pmfoc1a(:)   = 0.0_wp
[3124]1795
[3899]1796    IF ( init_aerosol_type == 1 )  THEN
[2505]1797!
[3864]1798!--    Read input profiles from PIDS_DYNAMIC_SALSA
[3483]1799#if defined( __netcdf )
[3864]1800!
1801!--    Location-dependent size distributions and compositions.
1802       INQUIRE( FILE = TRIM( input_file_dynamic ) //  TRIM( coupling_char ), EXIST = netcdf_extend )
[3124]1803       IF ( netcdf_extend )  THEN
1804!
[3864]1805!--       Open file in read-only mode
[3899]1806          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
[3124]1807!
[3864]1808!--       Inquire dimensions:
1809          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' )
1810          IF ( pr_nz /= nz )  THEN
1811             WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//&
1812                                        'the number of numeric grid points.'
1813             CALL message( 'aerosol_init', 'PA0601', 1, 2, 0, 6, 0 )
1814          ENDIF
1815          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_ncc, 'composition_index' )
1816!
1817!--       Allocate memory
[4012]1818          ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                              &
[3864]1819                    pr_mass_fracs_b(nzb:nzt+1,pr_ncc) )
1820          pr_mass_fracs_a = 0.0_wp
1821          pr_mass_fracs_b = 0.0_wp
1822!
1823!--       Read vertical levels
1824          CALL get_variable( id_dyn, 'z', pr_z )
1825!
1826!--       Read name and index of chemical components
1827          CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc )
1828          DO  ic = 1, pr_ncc
1829             SELECT CASE ( TRIM( cc_name(ic) ) )
1830                CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' )
1831                   cc_input_to_model(1) = ic
1832                CASE ( 'OC', 'oc' )
1833                   cc_input_to_model(2) = ic
1834                CASE ( 'BC', 'bc' )
1835                   cc_input_to_model(3) = ic
1836                CASE ( 'DU', 'du' )
1837                   cc_input_to_model(4) = ic
1838                CASE ( 'SS', 'ss' )
1839                   cc_input_to_model(5) = ic
1840                CASE ( 'HNO3', 'hno3', 'NO', 'no' )
1841                   cc_input_to_model(6) = ic
1842                CASE ( 'NH3', 'nh3', 'NH', 'nh' )
1843                   cc_input_to_model(7) = ic
1844             END SELECT
1845          ENDDO
1846
1847          IF ( SUM( cc_input_to_model ) == 0 )  THEN
1848             message_string = 'None of the aerosol chemical components in ' // TRIM(               &
1849                              input_file_dynamic ) // ' correspond to ones applied in SALSA.'
1850             CALL message( 'salsa_mod: aerosol_init', 'PA0602', 2, 2, 0, 6, 0 )
1851          ENDIF
1852!
1853!--       Vertical profiles of mass fractions of different chemical components:
1854          CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a,              &
1855                             0, pr_ncc-1, 0, pr_nz-1 )
1856          CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_b', pr_mass_fracs_b,              &
1857                             0, pr_ncc-1, 0, pr_nz-1  )
1858!
1859!--       Match the input data with the chemical composition applied in the model
1860          DO  ic = 1, maxspec
1861             ss = cc_input_to_model(ic)
1862             IF ( ss == 0 )  CYCLE
1863             pmf2a(nzb+1:nzt+1,ic) = pr_mass_fracs_a(nzb:nzt,ss)
1864             pmf2b(nzb+1:nzt+1,ic) = pr_mass_fracs_b(nzb:nzt,ss)
1865          ENDDO
1866!
1867!--       Aerosol concentrations: lod=1 (total PM) or lod=2 (sectional number size distribution)
1868          CALL get_attribute( id_dyn, 'lod', lod_aero, .FALSE., 'init_atmosphere_aerosol' )
[3899]1869          IF ( lod_aero /= 1 )  THEN
1870             message_string = 'Currently only lod=1 accepted for init_atmosphere_aerosol'
[3864]1871             CALL message( 'salsa_mod: aerosol_init', 'PA0603', 2, 2, 0, 6, 0 )
1872          ELSE
1873!
1874!--          Bin mean diameters in the input file
1875             CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nbins, 'Dmid')
1876             IF ( pr_nbins /= nbins_aerosol )  THEN
1877                message_string = 'Number of size bins in init_atmosphere_aerosol does not match '  &
1878                                 // 'with that applied in the model'
1879                CALL message( 'salsa_mod: aerosol_init', 'PA0604', 2, 2, 0, 6, 0 )
[3124]1880             ENDIF
[3864]1881
1882             ALLOCATE( pr_dmid(pr_nbins) )
1883             pr_dmid    = 0.0_wp
1884
1885             CALL get_variable( id_dyn, 'Dmid', pr_dmid )
1886!
1887!--          Check whether the sectional representation conform to the one
1888!--          applied in the model
1889             IF ( ANY( ABS( ( aero(1:nbins_aerosol)%dmid - pr_dmid ) /                             &
1890                              aero(1:nbins_aerosol)%dmid )  > 0.1_wp )  ) THEN
[3899]1891                message_string = 'Mean diameters of the aerosol size bins in ' // TRIM(            &
1892                                 input_file_dynamic ) // ' do not match with the sectional '//     &
[3864]1893                                 'representation of the model.'
1894                CALL message( 'salsa_mod: aerosol_init', 'PA0605', 2, 2, 0, 6, 0 )
[3124]1895             ENDIF
[3169]1896!
[3864]1897!--          Inital aerosol concentrations
1898             CALL get_variable( id_dyn, 'init_atmosphere_aerosol', pndist(nzb+1:nzt,:),            &
1899                                0, pr_nbins-1, 0, pr_nz-1 )
1900          ENDIF
1901!
1902!--       Set bottom and top boundary condition (Neumann)
1903          pmf2a(nzb,:)    = pmf2a(nzb+1,:)
1904          pmf2a(nzt+1,:)  = pmf2a(nzt,:)
1905          pmf2b(nzb,:)    = pmf2b(nzb+1,:)
1906          pmf2b(nzt+1,:)  = pmf2b(nzt,:)
1907          pndist(nzb,:)   = pndist(nzb+1,:)
1908          pndist(nzt+1,:) = pndist(nzt,:)
1909
1910          IF ( index_so4 < 0 )  THEN
1911             pmf2a(:,1) = 0.0_wp
1912             pmf2b(:,1) = 0.0_wp
1913          ENDIF
1914          IF ( index_oc < 0 )  THEN
1915             pmf2a(:,2) = 0.0_wp
1916             pmf2b(:,2) = 0.0_wp
1917          ENDIF
1918          IF ( index_bc < 0 )  THEN
1919             pmf2a(:,3) = 0.0_wp
1920             pmf2b(:,3) = 0.0_wp
1921          ENDIF
1922          IF ( index_du < 0 )  THEN
1923             pmf2a(:,4) = 0.0_wp
1924             pmf2b(:,4) = 0.0_wp
1925          ENDIF
1926          IF ( index_ss < 0 )  THEN
1927             pmf2a(:,5) = 0.0_wp
1928             pmf2b(:,5) = 0.0_wp
1929          ENDIF
1930          IF ( index_no < 0 )  THEN
1931             pmf2a(:,6) = 0.0_wp
1932             pmf2b(:,6) = 0.0_wp
1933          ENDIF
1934          IF ( index_nh < 0 )  THEN
1935             pmf2a(:,7) = 0.0_wp
1936             pmf2b(:,7) = 0.0_wp
1937          ENDIF
1938
1939          IF ( SUM( pmf2a ) < 0.00001_wp  .AND.  SUM( pmf2b ) < 0.00001_wp )  THEN
1940             message_string = 'Error in initialising mass fractions of chemical components. ' //   &
1941                              'Check that all chemical components are included in parameter file!'
1942             CALL message( 'salsa_mod: aerosol_init', 'PA0606', 2, 2, 0, 6, 0 ) 
1943          ENDIF
1944!
1945!--       Then normalise the mass fraction so that SUM = 1
1946          DO  k = nzb, nzt+1
1947             pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) )
1948             IF ( SUM( pmf2b(k,:) ) > 0.0_wp )  pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) )
1949          ENDDO
1950
1951          DEALLOCATE( pr_z, pr_mass_fracs_a, pr_mass_fracs_b )
1952
[3124]1953       ELSE
[3864]1954          message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) //  &
1955                           ' for SALSA missing!'
1956          CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 )
[4012]1957!
1958!--       Close input file
1959          CALL close_input_file( id_dyn )
[3864]1960       ENDIF   ! netcdf_extend
1961
1962#else
[3899]1963       message_string = 'init_aerosol_type = 1 but preprocessor directive __netcdf is not used '// &
1964                        'in compiling!'
[3864]1965       CALL message( 'salsa_mod: aerosol_init', 'PA0608', 1, 2, 0, 6, 0 )
1966
[3483]1967#endif
[3864]1968
[3899]1969    ELSEIF ( init_aerosol_type == 0 )  THEN
[3124]1970!
[3331]1971!--    Mass fractions for species in a and b-bins
[3864]1972       IF ( index_so4 > 0 )  THEN
1973          pmf2a(:,1) = mass_fracs_a(index_so4)
1974          pmf2b(:,1) = mass_fracs_b(index_so4)
[3331]1975       ENDIF
[3864]1976       IF ( index_oc > 0 )  THEN
1977          pmf2a(:,2) = mass_fracs_a(index_oc)
1978          pmf2b(:,2) = mass_fracs_b(index_oc)
[3331]1979       ENDIF
[3864]1980       IF ( index_bc > 0 )  THEN
1981          pmf2a(:,3) = mass_fracs_a(index_bc)
1982          pmf2b(:,3) = mass_fracs_b(index_bc)
[3331]1983       ENDIF
[3864]1984       IF ( index_du > 0 )  THEN
1985          pmf2a(:,4) = mass_fracs_a(index_du)
1986          pmf2b(:,4) = mass_fracs_b(index_du)
[3331]1987       ENDIF
[3864]1988       IF ( index_ss > 0 )  THEN
1989          pmf2a(:,5) = mass_fracs_a(index_ss)
1990          pmf2b(:,5) = mass_fracs_b(index_ss)
[3331]1991       ENDIF
[3864]1992       IF ( index_no > 0 )  THEN
1993          pmf2a(:,6) = mass_fracs_a(index_no)
1994          pmf2b(:,6) = mass_fracs_b(index_no)
[3331]1995       ENDIF
[3864]1996       IF ( index_nh > 0 )  THEN
1997          pmf2a(:,7) = mass_fracs_a(index_nh)
1998          pmf2b(:,7) = mass_fracs_b(index_nh)
1999       ENDIF
2000       DO  k = nzb, nzt+1
2001          pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) )
2002          IF ( SUM( pmf2b(k,:) ) > 0.0_wp ) pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) )
[3412]2003       ENDDO
[3864]2004
[3331]2005       CALL size_distribution( n_lognorm, dpg, sigmag, nsect )
2006!
2007!--    Normalize by the given total number concentration
[3864]2008       nsect = nsect * SUM( n_lognorm ) / SUM( nsect )
2009       DO  ib = start_subrange_1a, end_subrange_2b
2010          pndist(:,ib) = nsect(ib)
[3331]2011       ENDDO
2012    ENDIF
[3864]2013
[3899]2014    IF ( init_gases_type == 1 )  THEN
[3331]2015!
[3864]2016!--    Read input profiles from PIDS_CHEM
[3483]2017#if defined( __netcdf )
[3864]2018!
2019!--    Location-dependent size distributions and compositions.
2020       INQUIRE( FILE = TRIM( input_file_dynamic ) //  TRIM( coupling_char ), EXIST = netcdf_extend )
[3124]2021       IF ( netcdf_extend  .AND.  .NOT. salsa_gases_from_chem )  THEN
2022!
[3864]2023!--       Open file in read-only mode
[3899]2024          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
[3124]2025!
[3864]2026!--       Inquire dimensions:
2027          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' )
2028          IF ( pr_nz /= nz )  THEN
2029             WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//&
2030                                        'the number of numeric grid points.'
2031             CALL message( 'aerosol_init', 'PA0609', 1, 2, 0, 6, 0 )
2032          ENDIF
2033!
2034!--       Read vertical profiles of gases:
2035          CALL get_variable( id_dyn, 'init_atmosphere_h2so4', salsa_gas(1)%init(nzb+1:nzt) )
2036          CALL get_variable( id_dyn, 'init_atmosphere_hno3',  salsa_gas(2)%init(nzb+1:nzt) )
2037          CALL get_variable( id_dyn, 'init_atmosphere_nh3',   salsa_gas(3)%init(nzb+1:nzt) )
2038          CALL get_variable( id_dyn, 'init_atmosphere_ocnv',  salsa_gas(4)%init(nzb+1:nzt) )
2039          CALL get_variable( id_dyn, 'init_atmosphere_ocsv',  salsa_gas(5)%init(nzb+1:nzt) )
2040!
2041!--       Set Neumann top and surface boundary condition for initial + initialise concentrations
2042          DO  ig = 1, ngases_salsa
2043             salsa_gas(ig)%init(nzb)   =  salsa_gas(ig)%init(nzb+1)
2044             salsa_gas(ig)%init(nzt+1) =  salsa_gas(ig)%init(nzt)
[4012]2045             IF ( .NOT. read_restart_data_salsa )  THEN
2046                DO  k = nzb, nzt+1
2047                   salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
2048                ENDDO
2049             ENDIF
[3124]2050          ENDDO
[3864]2051
[3308]2052       ELSEIF ( .NOT. netcdf_extend  .AND.  .NOT.  salsa_gases_from_chem )  THEN
[3864]2053          message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) //  &
2054                           ' for SALSA missing!'
2055          CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 )
[4012]2056!
2057!--       Close input file
2058          CALL close_input_file( id_dyn )
[3864]2059       ENDIF   ! netcdf_extend
2060#else
[3899]2061       message_string = 'init_gases_type = 1 but preprocessor directive __netcdf is not used in '//&
2062                        'compiling!'
[3864]2063       CALL message( 'salsa_mod: aerosol_init', 'PA0611', 1, 2, 0, 6, 0 )
2064
[3483]2065#endif
2066
[2505]2067    ENDIF
[3864]2068!
2069!-- Both SO4 and OC are included, so use the given mass fractions
2070    IF ( index_oc > 0  .AND.  index_so4 > 0 )  THEN
2071       pmfoc1a(:) = pmf2a(:,2) / ( pmf2a(:,2) + pmf2a(:,1) )  ! Normalize
2072!
2073!-- Pure organic carbon
2074    ELSEIF ( index_oc > 0 )  THEN
2075       pmfoc1a(:) = 1.0_wp
2076!
2077!-- Pure SO4
2078    ELSEIF ( index_so4 > 0 )  THEN
2079       pmfoc1a(:) = 0.0_wp
[3169]2080
[3124]2081    ELSE
2082       message_string = 'Either OC or SO4 must be active for aerosol region 1a!'
[3864]2083       CALL message( 'salsa_mod: aerosol_init', 'PA0612', 1, 2, 0, 6, 0 )
2084    ENDIF
2085
[2505]2086!
2087!-- Initialize concentrations
[3864]2088    DO  i = nxlg, nxrg
[2505]2089       DO  j = nysg, nyng
[2777]2090          DO  k = nzb, nzt+1
[2505]2091!
[3864]2092!--          Predetermine flag to mask topography
[2505]2093             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[3864]2094!
[2505]2095!--          a) Number concentrations
[3864]2096!--          Region 1:
2097             DO  ib = start_subrange_1a, end_subrange_1a
[4012]2098                IF ( .NOT. read_restart_data_salsa )  THEN
2099                   aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
2100                ENDIF
[2744]2101                IF ( prunmode == 1 )  THEN
[3864]2102                   aerosol_number(ib)%init = pndist(:,ib)
[2744]2103                ENDIF
[2505]2104             ENDDO
[3864]2105!
2106!--          Region 2:
[2505]2107             IF ( nreg > 1 )  THEN
[3864]2108                DO  ib = start_subrange_2a, end_subrange_2a
[4012]2109                   IF ( .NOT. read_restart_data_salsa )  THEN
2110                      aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
2111                   ENDIF
[2744]2112                   IF ( prunmode == 1 )  THEN
[3864]2113                      aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib)
[2744]2114                   ENDIF
[2505]2115                ENDDO
2116                IF ( .NOT. no_insoluble )  THEN
[3864]2117                   DO  ib = start_subrange_2b, end_subrange_2b
2118                      IF ( pnf2a(k) < 1.0_wp )  THEN
[4012]2119                         IF ( .NOT. read_restart_data_salsa )  THEN
2120                            aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *    &
2121                                                             pndist(k,ib) * flag
2122                         ENDIF
[2744]2123                         IF ( prunmode == 1 )  THEN
[3864]2124                            aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp - nf2a ) * pndist(:,ib)
[2744]2125                         ENDIF
[2505]2126                      ENDIF
2127                   ENDDO
2128                ENDIF
2129             ENDIF
2130!
2131!--          b) Aerosol mass concentrations
[2744]2132!--             bin subrange 1: done here separately due to the SO4/OC convention
[3864]2133!
[2814]2134!--          SO4:
[3864]2135             IF ( index_so4 > 0 )  THEN
2136                ss = ( index_so4 - 1 ) * nbins_aerosol + start_subrange_1a !< start
2137                ee = ( index_so4 - 1 ) * nbins_aerosol + end_subrange_1a !< end
2138                ib = start_subrange_1a
2139                DO  ic = ss, ee
[4012]2140                   IF ( .NOT. read_restart_data_salsa )  THEN
2141                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) *          &
2142                                                     pndist(k,ib) * core(ib) * arhoh2so4 * flag
2143                   ENDIF
[2744]2144                   IF ( prunmode == 1 )  THEN
[3864]2145                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib) &
2146                                                 * core(ib) * arhoh2so4
[2744]2147                   ENDIF
[3864]2148                   ib = ib+1
[2505]2149                ENDDO
2150             ENDIF
[3864]2151!
[2814]2152!--          OC:
[3864]2153             IF ( index_oc > 0 ) THEN
2154                ss = ( index_oc - 1 ) * nbins_aerosol + start_subrange_1a !< start
2155                ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a !< end
2156                ib = start_subrange_1a
[4012]2157                DO  ic = ss, ee
2158                   IF ( .NOT. read_restart_data_salsa )  THEN
2159                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *    &
2160                                                     core(ib) * arhooc * flag
2161                   ENDIF
[2744]2162                   IF ( prunmode == 1 )  THEN
[3864]2163                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *        &
2164                                                 core(ib) * arhooc
[2744]2165                   ENDIF
[3864]2166                   ib = ib+1
[2505]2167                ENDDO
2168             ENDIF
2169          ENDDO !< k
[3864]2170
2171          prunmode = 3  ! Init only once
2172
[2505]2173       ENDDO !< j
2174    ENDDO !< i
[3864]2175
[2505]2176!
2177!-- c) Aerosol mass concentrations
[2744]2178!--    bin subrange 2:
[2505]2179    IF ( nreg > 1 ) THEN
[3864]2180
2181       IF ( index_so4 > 0 ) THEN
2182          CALL set_aero_mass( index_so4, pmf2a(:,1), pmf2b(:,1), pnf2a, pndist, core, arhoh2so4 )
[2505]2183       ENDIF
[3864]2184       IF ( index_oc > 0 ) THEN
2185          CALL set_aero_mass( index_oc, pmf2a(:,2), pmf2b(:,2), pnf2a, pndist, core, arhooc )
[2505]2186       ENDIF
[3864]2187       IF ( index_bc > 0 ) THEN
2188          CALL set_aero_mass( index_bc, pmf2a(:,3), pmf2b(:,3), pnf2a, pndist, core, arhobc )
[2505]2189       ENDIF
[3864]2190       IF ( index_du > 0 ) THEN
2191          CALL set_aero_mass( index_du, pmf2a(:,4), pmf2b(:,4), pnf2a, pndist, core, arhodu )
[2505]2192       ENDIF
[3864]2193       IF ( index_ss > 0 ) THEN
2194          CALL set_aero_mass( index_ss, pmf2a(:,5), pmf2b(:,5), pnf2a, pndist, core, arhoss )
[2505]2195       ENDIF
[3864]2196       IF ( index_no > 0 ) THEN
2197          CALL set_aero_mass( index_no, pmf2a(:,6), pmf2b(:,6), pnf2a, pndist, core, arhohno3 )
[2505]2198       ENDIF
[3864]2199       IF ( index_nh > 0 ) THEN
2200          CALL set_aero_mass( index_nh, pmf2a(:,7), pmf2b(:,7), pnf2a, pndist, core, arhonh3 )
[2505]2201       ENDIF
2202
2203    ENDIF
[3864]2204
[2505]2205 END SUBROUTINE aerosol_init
[3864]2206
[3308]2207!------------------------------------------------------------------------------!
2208! Description:
2209! ------------
2210!> Create a lognormal size distribution and discretise to a sectional
2211!> representation.
2212!------------------------------------------------------------------------------!
2213 SUBROUTINE size_distribution( in_ntot, in_dpg, in_sigma, psd_sect )
[3864]2214
[3308]2215    IMPLICIT NONE
[3864]2216
2217    INTEGER(iwp) ::  ib         !< running index: bin
2218    INTEGER(iwp) ::  iteration  !< running index: iteration
2219
2220    REAL(wp) ::  d1         !< particle diameter (m, dummy)
2221    REAL(wp) ::  d2         !< particle diameter (m, dummy)
2222    REAL(wp) ::  delta_d    !< (d2-d1)/10
2223    REAL(wp) ::  deltadp    !< bin width
2224    REAL(wp) ::  dmidi      !< ( d1 + d2 ) / 2
2225
2226    REAL(wp), DIMENSION(:), INTENT(in) ::  in_dpg    !< geometric mean diameter (m)
2227    REAL(wp), DIMENSION(:), INTENT(in) ::  in_ntot   !< number conc. (#/m3)
[3308]2228    REAL(wp), DIMENSION(:), INTENT(in) ::  in_sigma  !< standard deviation
[3864]2229
2230    REAL(wp), DIMENSION(:), INTENT(inout) ::  psd_sect  !< sectional size distribution
2231
2232    DO  ib = start_subrange_1a, end_subrange_2b
2233       psd_sect(ib) = 0.0_wp
2234!
[3308]2235!--    Particle diameter at the low limit (largest in the bin) (m)
[3864]2236       d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp
2237!
[3308]2238!--    Particle diameter at the high limit (smallest in the bin) (m)
[3864]2239       d2 = ( aero(ib)%vhilim / api6 )**0.33333333_wp
2240!
[3308]2241!--    Span of particle diameter in a bin (m)
[3864]2242       delta_d = 0.1_wp * ( d2 - d1 )
2243!
2244!--    Iterate:
2245       DO  iteration = 1, 10
2246          d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp + ( ib - 1) * delta_d
[3308]2247          d2 = d1 + delta_d
[3864]2248          dmidi = 0.5_wp * ( d1 + d2 )
[3308]2249          deltadp = LOG10( d2 / d1 )
[3864]2250!
[3308]2251!--       Size distribution
2252!--       in_ntot = total number, total area, or total volume concentration
2253!--       in_dpg = geometric-mean number, area, or volume diameter
2254!--       n(k) = number, area, or volume concentration in a bin
[3864]2255          psd_sect(ib) = psd_sect(ib) + SUM( in_ntot * deltadp / ( SQRT( 2.0_wp * pi ) *           &
2256                        LOG10( in_sigma ) ) * EXP( -LOG10( dmidi / in_dpg )**2.0_wp /              &
2257                        ( 2.0_wp * LOG10( in_sigma ) ** 2.0_wp ) ) )
2258
[3308]2259       ENDDO
2260    ENDDO
[3864]2261
[3308]2262 END SUBROUTINE size_distribution
[2505]2263
2264!------------------------------------------------------------------------------!
2265! Description:
2266! ------------
2267!> Sets the mass concentrations to aerosol arrays in 2a and 2b.
2268!>
2269!> Tomi Raatikainen, FMI, 29.2.2016
2270!------------------------------------------------------------------------------!
[3864]2271 SUBROUTINE set_aero_mass( ispec, pmf2a, pmf2b, pnf2a, pndist, pcore, prho )
2272
[2505]2273    IMPLICIT NONE
2274
[3864]2275    INTEGER(iwp) ::  ee        !< index: end
2276    INTEGER(iwp) ::  i         !< loop index
2277    INTEGER(iwp) ::  ib        !< loop index
2278    INTEGER(iwp) ::  ic        !< loop index
2279    INTEGER(iwp) ::  j         !< loop index
2280    INTEGER(iwp) ::  k         !< loop index
2281    INTEGER(iwp) ::  prunmode  !< 1 = initialise
2282    INTEGER(iwp) ::  ss        !< index: start
2283
[2505]2284    INTEGER(iwp), INTENT(in) :: ispec  !< Aerosol species index
[3864]2285
2286    REAL(wp) ::  flag   !< flag to mask topography grid points
2287
[2505]2288    REAL(wp), INTENT(in) ::  prho !< Aerosol density
[3864]2289
2290    REAL(wp), DIMENSION(nbins_aerosol), INTENT(in) ::  pcore !< Aerosol bin mid core volume
2291    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pnf2a !< Number fraction for 2a
2292    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pmf2a !< Mass distributions for a
2293    REAL(wp), DIMENSION(0:nz+1), INTENT(in)        ::  pmf2b !< and b bins
2294
2295    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol), INTENT(in) ::  pndist !< Aerosol size distribution
2296
[2744]2297    prunmode = 1
[3864]2298
2299    DO i = nxlg, nxrg
[2505]2300       DO j = nysg, nyng
[3864]2301          DO k = nzb, nzt+1
[2505]2302!
2303!--          Predetermine flag to mask topography
2304             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
[3864]2305!
[2505]2306!--          Regime 2a:
[3864]2307             ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2a
2308             ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2a
2309             ib = start_subrange_2a
2310             DO ic = ss, ee
[4012]2311                IF ( .NOT. read_restart_data_salsa )  THEN
2312                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib)&
2313                                                  * pcore(ib) * prho * flag
2314                ENDIF
[2744]2315                IF ( prunmode == 1 )  THEN
[3864]2316                   aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) *  &
2317                                              pcore(ib) * prho
[2744]2318                ENDIF
[3864]2319                ib = ib + 1
[2505]2320             ENDDO
[3864]2321!
[2505]2322!--          Regime 2b:
2323             IF ( .NOT. no_insoluble )  THEN
[3864]2324                ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2b
2325                ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2b
2326                ib = start_subrange_2a
2327                DO ic = ss, ee
[4012]2328                   IF ( .NOT. read_restart_data_salsa )  THEN
2329                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k))&
2330                                                     * pndist(k,ib) * pcore(ib) * prho * flag
2331                   ENDIF
[2744]2332                   IF ( prunmode == 1 )  THEN
[3864]2333                      aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) * &
2334                                                 pndist(k,ib) * pcore(ib) * prho
[2744]2335                   ENDIF
[3864]2336                   ib = ib + 1
2337                ENDDO  ! c
2338
[2505]2339             ENDIF
[3864]2340          ENDDO   ! k
2341
2342          prunmode = 3  ! Init only once
2343
2344       ENDDO   ! j
2345    ENDDO   ! i
2346
2347 END SUBROUTINE set_aero_mass
2348
2349!------------------------------------------------------------------------------!
2350! Description:
2351! ------------
2352!> Initialise the matching between surface types in LSM and deposition models.
2353!> Do the matching based on Zhang et al. (2001). Atmos. Environ. 35, 549-560
2354!> (here referred as Z01).
2355!------------------------------------------------------------------------------!
2356 SUBROUTINE init_deposition
2357
2358    USE surface_mod,                                                                               &
[3956]2359        ONLY:  surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[3864]2360
2361    IMPLICIT NONE
2362
2363    INTEGER(iwp) ::  l  !< loop index for vertical surfaces
2364
[3956]2365    LOGICAL :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM surfaces)
[3864]2366
[3956]2367    IF ( depo_pcm_par == 'zhang2001' )  THEN
2368       depo_pcm_par_num = 1
2369    ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
2370       depo_pcm_par_num = 2
2371    ENDIF
[3864]2372
[3956]2373    IF ( depo_surf_par == 'zhang2001' )  THEN
2374       depo_surf_par_num = 1
2375    ELSEIF ( depo_surf_par == 'petroff2010' )  THEN
2376       depo_surf_par_num = 2
2377    ENDIF
2378!
2379!-- LSM: Pavement, vegetation and water
2380    IF ( nldepo_surf  .AND.  land_surface )  THEN
2381       match_lsm = .TRUE.
2382       ALLOCATE( lsm_to_depo_h%match_lupg(1:surf_lsm_h%ns),                                         &
2383                 lsm_to_depo_h%match_luvw(1:surf_lsm_h%ns),                                         &
2384                 lsm_to_depo_h%match_luww(1:surf_lsm_h%ns) )
2385       lsm_to_depo_h%match_lupg = 0
2386       lsm_to_depo_h%match_luvw = 0
2387       lsm_to_depo_h%match_luww = 0
2388       CALL match_sm_zhang( surf_lsm_h, lsm_to_depo_h%match_lupg, lsm_to_depo_h%match_luvw,        &
2389                            lsm_to_depo_h%match_luww, match_lsm )
[3864]2390       DO  l = 0, 3
[3956]2391          ALLOCATE( lsm_to_depo_v(l)%match_lupg(1:surf_lsm_v(l)%ns),                               &
2392                    lsm_to_depo_v(l)%match_luvw(1:surf_lsm_v(l)%ns),                               &
2393                    lsm_to_depo_v(l)%match_luww(1:surf_lsm_v(l)%ns) )
2394          lsm_to_depo_v(l)%match_lupg = 0
2395          lsm_to_depo_v(l)%match_luvw = 0
2396          lsm_to_depo_v(l)%match_luww = 0
2397          CALL match_sm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match_lupg,                         &
2398                               lsm_to_depo_v(l)%match_luvw, lsm_to_depo_v(l)%match_luww, match_lsm )
[2505]2399       ENDDO
[3864]2400    ENDIF
[3956]2401!
2402!-- USM: Green roofs/walls, wall surfaces and windows
2403    IF ( nldepo_surf  .AND.  urban_surface )  THEN
2404       match_lsm = .FALSE.
2405       ALLOCATE( usm_to_depo_h%match_lupg(1:surf_usm_h%ns),                                        &
2406                 usm_to_depo_h%match_luvw(1:surf_usm_h%ns),                                        &
2407                 usm_to_depo_h%match_luww(1:surf_usm_h%ns) )
2408       usm_to_depo_h%match_lupg = 0
2409       usm_to_depo_h%match_luvw = 0
2410       usm_to_depo_h%match_luww = 0
2411       CALL match_sm_zhang( surf_usm_h, usm_to_depo_h%match_lupg, usm_to_depo_h%match_luvw,        &
2412                            usm_to_depo_h%match_luww, match_lsm )
2413       DO  l = 0, 3
2414          ALLOCATE( usm_to_depo_v(l)%match_lupg(1:surf_usm_v(l)%ns),                               &
2415                    usm_to_depo_v(l)%match_luvw(1:surf_usm_v(l)%ns),                               &
2416                    usm_to_depo_v(l)%match_luww(1:surf_usm_v(l)%ns) )
2417          usm_to_depo_v(l)%match_lupg = 0
2418          usm_to_depo_v(l)%match_luvw = 0
2419          usm_to_depo_v(l)%match_luww = 0
2420          CALL match_sm_zhang( surf_usm_v(l), usm_to_depo_v(l)%match_lupg,                         &
2421                               usm_to_depo_v(l)%match_luvw, usm_to_depo_v(l)%match_luww, match_lsm )
2422       ENDDO
2423    ENDIF
[3864]2424
2425    IF ( nldepo_pcm )  THEN
2426       SELECT CASE ( depo_pcm_type )
2427          CASE ( 'evergreen_needleleaf' )
2428             depo_pcm_type_num = 1
2429          CASE ( 'evergreen_broadleaf' )
2430             depo_pcm_type_num = 2
2431          CASE ( 'deciduous_needleleaf' )
2432             depo_pcm_type_num = 3
2433          CASE ( 'deciduous_broadleaf' )
2434             depo_pcm_type_num = 4
2435          CASE DEFAULT
2436             message_string = 'depo_pcm_type not set correctly.'
2437             CALL message( 'salsa_mod: init_deposition', 'PA0613', 1, 2, 0, 6, 0 )
2438       END SELECT
2439    ENDIF
2440
2441 END SUBROUTINE init_deposition
2442
2443!------------------------------------------------------------------------------!
2444! Description:
2445! ------------
2446!> Match the surface types in PALM and Zhang et al. 2001 deposition module
2447!------------------------------------------------------------------------------!
[3956]2448 SUBROUTINE match_sm_zhang( surf, match_pav_green, match_veg_wall, match_wat_win, match_lsm )
[3864]2449
2450    USE surface_mod,                                                           &
2451        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
2452
2453    IMPLICIT NONE
2454
[3956]2455    INTEGER(iwp) ::  m              !< index for surface elements
2456    INTEGER(iwp) ::  pav_type_palm  !< pavement / green wall type in PALM
2457    INTEGER(iwp) ::  veg_type_palm  !< vegetation / wall type in PALM
2458    INTEGER(iwp) ::  wat_type_palm  !< water / window type in PALM
[3864]2459
[3956]2460    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_pav_green  !<  matching pavement/green walls
2461    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_veg_wall   !<  matching vegetation/walls
2462    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_wat_win    !<  matching water/windows
2463
2464    LOGICAL, INTENT(in) :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM)
2465
[3864]2466    TYPE(surf_type), INTENT(in) :: surf  !< respective surface type
2467
2468    DO  m = 1, surf%ns
[3956]2469       IF ( match_lsm )  THEN
2470!
2471!--       Vegetation (LSM):
2472          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
2473             veg_type_palm = surf%vegetation_type(m)
2474             SELECT CASE ( veg_type_palm )
2475                CASE ( 0 )
2476                   message_string = 'No vegetation type defined.'
2477                   CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
2478                CASE ( 1 )  ! bare soil
2479                   match_veg_wall(m) = 6  ! grass in Z01
2480                CASE ( 2 )  ! crops, mixed farming
2481                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
2482                CASE ( 3 )  ! short grass
2483                   match_veg_wall(m) = 6  ! grass in Z01
2484                CASE ( 4 )  ! evergreen needleleaf trees
2485                    match_veg_wall(m) = 1  ! evergreen needleleaf trees in Z01
2486                CASE ( 5 )  ! deciduous needleleaf trees
2487                   match_veg_wall(m) = 3  ! deciduous needleleaf trees in Z01
2488                CASE ( 6 )  ! evergreen broadleaf trees
2489                   match_veg_wall(m) = 2  ! evergreen broadleaf trees in Z01
2490                CASE ( 7 )  ! deciduous broadleaf trees
2491                   match_veg_wall(m) = 4  ! deciduous broadleaf trees in Z01
2492                CASE ( 8 )  ! tall grass
2493                   match_veg_wall(m) = 6  ! grass in Z01
2494                CASE ( 9 )  ! desert
2495                   match_veg_wall(m) = 8  ! desert in Z01
2496                CASE ( 10 )  ! tundra
2497                   match_veg_wall(m) = 9  ! tundra in Z01
2498                CASE ( 11 )  ! irrigated crops
2499                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
2500                CASE ( 12 )  ! semidesert
2501                   match_veg_wall(m) = 8  ! desert in Z01
2502                CASE ( 13 )  ! ice caps and glaciers
2503                   match_veg_wall(m) = 12  ! ice cap and glacier in Z01
2504                CASE ( 14 )  ! bogs and marshes
2505                   match_veg_wall(m) = 11  ! wetland with plants in Z01
2506                CASE ( 15 )  ! evergreen shrubs
2507                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2508                CASE ( 16 )  ! deciduous shrubs
2509                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2510                CASE ( 17 )  ! mixed forest/woodland
2511                   match_veg_wall(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
2512                CASE ( 18 )  ! interrupted forest
2513                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
2514             END SELECT
[3864]2515          ENDIF
[3956]2516!
2517!--       Pavement (LSM):
2518          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
2519             pav_type_palm = surf%pavement_type(m)
2520             IF ( pav_type_palm == 0 )  THEN  ! error
2521                message_string = 'No pavement type defined.'
2522                CALL message( 'salsa_mod: match_sm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
2523             ELSE
2524                match_pav_green(m) = 15  ! urban in Z01
2525             ENDIF
[3864]2526          ENDIF
[3956]2527!
2528!--       Water (LSM):
2529          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
2530             wat_type_palm = surf%water_type(m)
2531             IF ( wat_type_palm == 0 )  THEN  ! error
2532                message_string = 'No water type defined.'
2533                CALL message( 'salsa_mod: match_sm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
2534             ELSEIF ( wat_type_palm == 3 )  THEN
2535                match_wat_win(m) = 14  ! ocean in Z01
2536             ELSEIF ( wat_type_palm == 1  .OR.  wat_type_palm == 2 .OR.  wat_type_palm == 4        &
2537                      .OR.  wat_type_palm == 5  )  THEN
2538                match_wat_win(m) = 13  ! inland water in Z01
2539             ENDIF
2540          ENDIF
2541       ELSE
2542!
2543!--       Wall surfaces (USM):
2544          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
2545             match_veg_wall(m) = 15  ! urban in Z01
2546          ENDIF
2547!
2548!--       Green walls and roofs (USM):
2549          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
2550             match_pav_green(m) =  6 ! (short) grass in Z01
2551          ENDIF
2552!
2553!--       Windows (USM):
2554          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
2555             match_wat_win(m) = 15  ! urban in Z01
2556          ENDIF
[3864]2557       ENDIF
2558
[2505]2559    ENDDO
2560
[3956]2561 END SUBROUTINE match_sm_zhang
[3864]2562
[2505]2563!------------------------------------------------------------------------------!
2564! Description:
2565! ------------
2566!> Swapping of timelevels
2567!------------------------------------------------------------------------------!
2568 SUBROUTINE salsa_swap_timelevel( mod_count )
2569
2570    IMPLICIT NONE
2571
[3864]2572    INTEGER(iwp) ::  ib   !<
2573    INTEGER(iwp) ::  ic   !<
2574    INTEGER(iwp) ::  icc  !<
2575    INTEGER(iwp) ::  ig   !<
2576
[2505]2577    INTEGER(iwp), INTENT(IN) ::  mod_count  !<
2578
[3864]2579    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
[3636]2580
[3864]2581       SELECT CASE ( mod_count )
[2505]2582
[3864]2583          CASE ( 0 )
[2505]2584
[3864]2585             DO  ib = 1, nbins_aerosol
2586                aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => nconc_1(:,:,:,ib)
2587                aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,ib)
2588
2589                DO  ic = 1, ncomponents_mass
2590                   icc = ( ic-1 ) * nbins_aerosol + ib
2591                   aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => mconc_1(:,:,:,icc)
2592                   aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,icc)
2593                ENDDO
[2505]2594             ENDDO
2595
[3864]2596             IF ( .NOT. salsa_gases_from_chem )  THEN
2597                DO  ig = 1, ngases_salsa
2598                   salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => gconc_1(:,:,:,ig)
2599                   salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,ig)
2600                ENDDO
2601             ENDIF
[2505]2602
[3864]2603          CASE ( 1 )
2604
2605             DO  ib = 1, nbins_aerosol
2606                aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => nconc_2(:,:,:,ib)
2607                aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,ib)
2608                DO  ic = 1, ncomponents_mass
2609                   icc = ( ic-1 ) * nbins_aerosol + ib
2610                   aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => mconc_2(:,:,:,icc)
2611                   aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,icc)
2612                ENDDO
[2505]2613             ENDDO
2614
[3864]2615             IF ( .NOT. salsa_gases_from_chem )  THEN
2616                DO  ig = 1, ngases_salsa
2617                   salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   => gconc_2(:,:,:,ig)
2618                   salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,ig)
2619                ENDDO
2620             ENDIF
[2505]2621
[3864]2622       END SELECT
2623
[3637]2624    ENDIF
2625
[2505]2626 END SUBROUTINE salsa_swap_timelevel
2627
2628
2629!------------------------------------------------------------------------------!
2630! Description:
2631! ------------
2632!> This routine reads the respective restart data.
2633!------------------------------------------------------------------------------!
[3864]2634 SUBROUTINE salsa_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,      &
2635                             nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
[2505]2636
[4117]2637    USE control_parameters,                                                                        &
2638        ONLY:  length, restart_string
2639
[2505]2640    IMPLICIT NONE
[3864]2641
2642    INTEGER(iwp) ::  ib              !<
2643    INTEGER(iwp) ::  ic              !<
2644    INTEGER(iwp) ::  ig              !<
2645    INTEGER(iwp) ::  k               !<
[3582]2646    INTEGER(iwp) ::  nxlc            !<
2647    INTEGER(iwp) ::  nxlf            !<
2648    INTEGER(iwp) ::  nxl_on_file     !<
2649    INTEGER(iwp) ::  nxrc            !<
2650    INTEGER(iwp) ::  nxrf            !<
2651    INTEGER(iwp) ::  nxr_on_file     !<
2652    INTEGER(iwp) ::  nync            !<
2653    INTEGER(iwp) ::  nynf            !<
2654    INTEGER(iwp) ::  nyn_on_file     !<
2655    INTEGER(iwp) ::  nysc            !<
2656    INTEGER(iwp) ::  nysf            !<
2657    INTEGER(iwp) ::  nys_on_file     !<
2658
[3864]2659    LOGICAL, INTENT(OUT)  ::  found  !<
[3582]2660
2661    REAL(wp), &
2662       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
[3864]2663
[3582]2664    found = .FALSE.
[3864]2665
[3005]2666    IF ( read_restart_data_salsa )  THEN
[3864]2667
[3582]2668       SELECT CASE ( restart_string(1:length) )
[3864]2669
[3582]2670          CASE ( 'aerosol_number' )
[3864]2671             DO  ib = 1, nbins_aerosol
[3582]2672                IF ( k == 1 )  READ ( 13 ) tmp_3d
[3864]2673                aerosol_number(ib)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
2674                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
[3582]2675                found = .TRUE.
[2505]2676             ENDDO
[3864]2677
[3582]2678          CASE ( 'aerosol_mass' )
[3864]2679             DO  ic = 1, ncomponents_mass * nbins_aerosol
[3582]2680                IF ( k == 1 )  READ ( 13 ) tmp_3d
[3864]2681                aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
2682                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
[3582]2683                found = .TRUE.
2684             ENDDO
[3864]2685
[3582]2686          CASE ( 'salsa_gas' )
[3864]2687             DO  ig = 1, ngases_salsa
[3582]2688                IF ( k == 1 )  READ ( 13 ) tmp_3d
[3864]2689                salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                    &
2690                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
[3582]2691                found = .TRUE.
2692             ENDDO
[3864]2693
[3582]2694          CASE DEFAULT
2695             found = .FALSE.
[3864]2696
[3582]2697       END SELECT
[2505]2698    ENDIF
2699
[3006]2700 END SUBROUTINE salsa_rrd_local
[2505]2701
2702!------------------------------------------------------------------------------!
2703! Description:
2704! ------------
2705!> This routine writes the respective restart data.
[2525]2706!> Note that the following input variables in PARIN have to be equal between
2707!> restart runs:
[2777]2708!>    listspec, nbin, nbin2, nf2a, ncc, mass_fracs_a, mass_fracs_b
[2505]2709!------------------------------------------------------------------------------!
[3006]2710 SUBROUTINE salsa_wrd_local
[2505]2711
[4117]2712    USE control_parameters,                                                                        &
2713        ONLY:  write_binary
2714
[2505]2715    IMPLICIT NONE
[3864]2716
2717    INTEGER(iwp) ::  ib   !<
2718    INTEGER(iwp) ::  ic   !<
2719    INTEGER(iwp) ::  ig  !<
2720
[3005]2721    IF ( write_binary  .AND.  write_binary_salsa )  THEN
[3864]2722
[3582]2723       CALL wrd_write_string( 'aerosol_number' )
[3864]2724       DO  ib = 1, nbins_aerosol
2725          WRITE ( 14 )  aerosol_number(ib)%conc
[3582]2726       ENDDO
[3864]2727
[3582]2728       CALL wrd_write_string( 'aerosol_mass' )
[3864]2729       DO  ic = 1, nbins_aerosol * ncomponents_mass
2730          WRITE ( 14 )  aerosol_mass(ic)%conc
[2525]2731       ENDDO
[3864]2732
[3582]2733       CALL wrd_write_string( 'salsa_gas' )
[3864]2734       DO  ig = 1, ngases_salsa
2735          WRITE ( 14 )  salsa_gas(ig)%conc
[2505]2736       ENDDO
[3864]2737
[2505]2738    ENDIF
2739
[3864]2740 END SUBROUTINE salsa_wrd_local
[2505]2741
2742!------------------------------------------------------------------------------!
2743! Description:
2744! ------------
2745!> Performs necessary unit and dimension conversion between the host model and
2746!> SALSA module, and calls the main SALSA routine.
2747!> Partially adobted form the original SALSA boxmodel version.
2748!> Now takes masses in as kg/kg from LES!! Converted to m3/m3 for SALSA
2749!> 05/2016 Juha: This routine is still pretty much in its original shape.
2750!>               It's dumb as a mule and twice as ugly, so implementation of
2751!>               an improved solution is necessary sooner or later.
2752!> Juha Tonttila, FMI, 2014
2753!> Jaakko Ahola, FMI, 2016
2754!> Only aerosol processes included, Mona Kurppa, UHel, 2017
2755!------------------------------------------------------------------------------!
[3169]2756 SUBROUTINE salsa_driver( i, j, prunmode )
[2505]2757
[3864]2758    USE arrays_3d,                                                                                 &
2759        ONLY: pt_p, q_p, u, v, w
2760
2761    USE plant_canopy_model_mod,                                                                    &
[3064]2762        ONLY: lad_s
[3864]2763
2764    USE surface_mod,                                                                               &
2765        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
2766
[2505]2767    IMPLICIT NONE
[3864]2768
2769    INTEGER(iwp) ::  endi    !< end index
2770    INTEGER(iwp) ::  ib      !< loop index
2771    INTEGER(iwp) ::  ic      !< loop index
2772    INTEGER(iwp) ::  ig      !< loop index
2773    INTEGER(iwp) ::  k_wall  !< vertical index of topography top
2774    INTEGER(iwp) ::  k       !< loop index
2775    INTEGER(iwp) ::  l       !< loop index
2776    INTEGER(iwp) ::  nc_h2o  !< index of H2O in the prtcl index table
2777    INTEGER(iwp) ::  ss      !< loop index
2778    INTEGER(iwp) ::  str     !< start index
2779    INTEGER(iwp) ::  vc      !< default index in prtcl
2780
2781    INTEGER(iwp), INTENT(in) ::  i         !< loop index
2782    INTEGER(iwp), INTENT(in) ::  j         !< loop index
2783    INTEGER(iwp), INTENT(in) ::  prunmode  !< 1: Initialization, 2: Spinup, 3: Regular runtime
2784
2785    REAL(wp) ::  cw_old  !< previous H2O mixing ratio
2786    REAL(wp) ::  flag    !< flag to mask topography grid points
2787    REAL(wp) ::  in_lad  !< leaf area density (m2/m3)
2788    REAL(wp) ::  in_rh   !< relative humidity
2789    REAL(wp) ::  zgso4   !< SO4
2790    REAL(wp) ::  zghno3  !< HNO3
2791    REAL(wp) ::  zgnh3   !< NH3
2792    REAL(wp) ::  zgocnv  !< non-volatile OC
2793    REAL(wp) ::  zgocsv  !< semi-volatile OC
2794
2795    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_adn  !< air density (kg/m3)
2796    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cs   !< H2O sat. vapour conc.
2797    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cw   !< H2O vapour concentration
2798    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_p    !< pressure (Pa)
2799    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_t    !< temperature (K)
2800    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_u    !< wind magnitude (m/s)
2801    REAL(wp), DIMENSION(nzb:nzt+1) ::  kvis    !< kinematic viscosity of air(m2/s)
2802    REAL(wp), DIMENSION(nzb:nzt+1) ::  ppm_to_nconc  !< Conversion factor from ppm to #/m3
2803
2804    REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) ::  schmidt_num  !< particle Schmidt number
2805    REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) ::  vd           !< particle fall seed (m/s)
2806
[3899]2807    TYPE(t_section), DIMENSION(nbins_aerosol) ::  lo_aero   !< additional variable for OpenMP
2808    TYPE(t_section), DIMENSION(nbins_aerosol) ::  aero_old  !< helper array
[3864]2809
[2505]2810    aero_old(:)%numc = 0.0_wp
[2754]2811    in_lad           = 0.0_wp
2812    in_u             = 0.0_wp
[3064]2813    kvis             = 0.0_wp
[3899]2814    lo_aero          = aero
[3864]2815    schmidt_num      = 0.0_wp
[3064]2816    vd               = 0.0_wp
[2754]2817    zgso4            = nclim
2818    zghno3           = nclim
2819    zgnh3            = nclim
2820    zgocnv           = nclim
2821    zgocsv           = nclim
[3864]2822!
[2505]2823!-- Aerosol number is always set, but mass can be uninitialized
[3864]2824    DO ib = 1, nbins_aerosol
[3899]2825       lo_aero(ib)%volc(:)  = 0.0_wp
[3864]2826       aero_old(ib)%volc(:) = 0.0_wp
2827    ENDDO
2828!
[2505]2829!-- Set the salsa runtime config (How to make this more efficient?)
2830    CALL set_salsa_runtime( prunmode )
[3864]2831!
[2505]2832!-- Calculate thermodynamic quantities needed in SALSA
[3864]2833    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 )
[2505]2834!
[3064]2835!-- Magnitude of wind: needed for deposition
[3169]2836    IF ( lsdepo )  THEN
[3864]2837       in_u(nzb+1:nzt) = SQRT( ( 0.5_wp * ( u(nzb+1:nzt,j,i) + u(nzb+1:nzt,j,i+1) ) )**2 +         &
2838                               ( 0.5_wp * ( v(nzb+1:nzt,j,i) + v(nzb+1:nzt,j+1,i) ) )**2 +         &
2839                               ( 0.5_wp * ( w(nzb:nzt-1,j,i) + w(nzb+1:nzt,j,  i) ) )**2 )
[3064]2840    ENDIF
2841!
[2754]2842!-- Calculate conversion factors for gas concentrations
[3864]2843    ppm_to_nconc(:) = for_ppm_to_nconc * in_p(:) / in_t(:)
[2754]2844!
[2505]2845!-- Determine topography-top index on scalar grid
[3864]2846    k_wall = k_topo_top(j,i)
2847
[3169]2848    DO k = nzb+1, nzt
[2505]2849!
2850!--    Predetermine flag to mask topography
[3169]2851       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2505]2852!
[3864]2853!--    Wind velocity for dry depositon on vegetation
2854       IF ( lsdepo_pcm  .AND.  plant_canopy )  THEN
2855          in_lad = lad_s( MAX( k-k_wall,0 ),j,i)
2856       ENDIF
[2814]2857!
[3018]2858!--    For initialization and spinup, limit the RH with the parameter rhlim
[2505]2859       IF ( prunmode < 3 ) THEN
[3169]2860          in_cw(k) = MIN( in_cw(k), in_cs(k) * rhlim )
[2505]2861       ELSE
[3169]2862          in_cw(k) = in_cw(k)
[2505]2863       ENDIF
[3169]2864       cw_old = in_cw(k) !* in_adn(k)
[3864]2865!
[2505]2866!--    Set volume concentrations:
2867!--    Sulphate (SO4) or sulphuric acid H2SO4
[3864]2868       IF ( index_so4 > 0 )  THEN
[2505]2869          vc = 1
[3864]2870          str = ( index_so4-1 ) * nbins_aerosol + 1    ! start index
2871          endi = index_so4 * nbins_aerosol             ! end index
2872          ic = 1
[2505]2873          DO ss = str, endi
[3899]2874             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2so4
[3864]2875             ic = ic+1
[2505]2876          ENDDO
[3899]2877          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
[2505]2878       ENDIF
[3864]2879!
[2505]2880!--    Organic carbon (OC) compounds
[3864]2881       IF ( index_oc > 0 )  THEN
[2505]2882          vc = 2
[3864]2883          str = ( index_oc-1 ) * nbins_aerosol + 1
2884          endi = index_oc * nbins_aerosol
2885          ic = 1
[2505]2886          DO ss = str, endi
[3899]2887             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhooc
[3864]2888             ic = ic+1
[2505]2889          ENDDO
[3899]2890          aero_old(1:nbins_aerosol)%volc(vc) = lo_aero(1:nbins_aerosol)%volc(vc)
[2505]2891       ENDIF
[3864]2892!
[2505]2893!--    Black carbon (BC)
[3864]2894       IF ( index_bc > 0 )  THEN
[2505]2895          vc = 3
[3864]2896          str = ( index_bc-1 ) * nbins_aerosol + 1 + end_subrange_1a
2897          endi = index_bc * nbins_aerosol
2898          ic = 1 + end_subrange_1a
[2505]2899          DO ss = str, endi
[3899]2900             lo_aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhobc
[3864]2901             ic = ic+1
2902          ENDDO