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

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

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

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