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

Last change on this file since 3787 was 3787, checked in by raasch, 4 years ago

unused variables removed

  • Property svn:keywords set to Id
File size: 458.2 KB
Line 
1!> @file salsa_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM-4U.
4!
5! PALM-4U is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2018-2018 University of Helsinki
18! Copyright 1997-2019 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: salsa_mod.f90 3787 2019-03-07 08:43:54Z raasch $
28! unused variables removed
29!
30! 3780 2019-03-05 11:19:45Z forkel
31! unused variable for file index removed from rrd-subroutines parameter list
32!
33! 3685 2019-01-21 01:02:11Z knoop
34! Some interface calls moved to module_interface + cleanup
35!
36! 3655 2019-01-07 16:51:22Z knoop
37! Implementation of the PALM module interface
38!
39! 3636 2018-12-19 13:48:34Z raasch
40! nopointer option removed
41!
42! 3630 2018-12-17 11:04:17Z knoop
43! - Moved the control parameter "salsa" from salsa_mod.f90 to control_parameters
44! - Updated salsa_rrd_local and salsa_wrd_local
45! - Add target attribute
46! - Revise initialization in case of restarts
47! - Revise masked data output
48!
49! 3582 2018-11-29 19:16:36Z suehring
50! missing comma separator inserted
51!
52! 3483 2018-11-02 14:19:26Z raasch
53! bugfix: directives added to allow compilation without netCDF
54!
55! 3481 2018-11-02 09:14:13Z raasch
56! temporary variable cc introduced to circumvent a possible Intel18 compiler bug
57! related to contiguous/non-contguous pointer/target attributes
58!
59! 3473 2018-10-30 20:50:15Z suehring
60! NetCDF input routine renamed
61!
62! 3467 2018-10-30 19:05:21Z suehring
63! Initial revision
64!
65! 3412 2018-10-24 07:25:57Z monakurppa
66!
67! Authors:
68! --------
69! @author Mona Kurppa (University of Helsinki)
70!
71!
72! Description:
73! ------------
74!> Sectional aerosol module for large scale applications SALSA
75!> (Kokkola et al., 2008, ACP 8, 2469-2483). Solves the aerosol number and mass
76!> concentration as well as chemical composition. Includes aerosol dynamic
77!> processes: nucleation, condensation/evaporation of vapours, coagulation and
78!> deposition on tree leaves, ground and roofs.
79!> Implementation is based on formulations implemented in UCLALES-SALSA except
80!> for deposition which is based on parametrisations by Zhang et al. (2001,
81!> Atmos. Environ. 35, 549-560) or Petroff&Zhang (2010, Geosci. Model Dev. 3,
82!> 753-769)
83!>
84!> @todo Implement turbulent inflow of aerosols in inflow_turbulence.
85!> @todo Deposition on subgrid scale vegetation
86!> @todo Deposition on vegetation calculated by default for deciduous broadleaf
87!>       trees
88!> @todo Revise masked data output. There is a potential bug in case of
89!>       terrain-following masked output, according to data_output_mask.
90!> @todo There are now improved interfaces for NetCDF data input which can be
91!>       used instead of get variable etc.
92!------------------------------------------------------------------------------!
93 MODULE salsa_mod
94
95    USE basic_constants_and_equations_mod,                                     &
96        ONLY:  c_p, g, p_0, pi, r_d
97 
98    USE chemistry_model_mod,                                                   &
99        ONLY:  chem_species, nspec, nvar, spc_names
100
101    USE chem_modules,                                                          &
102        ONLY:  call_chem_at_all_substeps, chem_gasphase_on
103
104    USE control_parameters
105
106    USE indices,                                                               &
107        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
108               nzb_s_inner, nz, nzt, wall_flags_0
109     
110    USE kinds
111   
112    USE pegrid
113   
114    USE salsa_util_mod
115
116    IMPLICIT NONE
117!
118!-- SALSA constants:
119!
120!-- Local constants:
121    INTEGER(iwp), PARAMETER ::  ngast   = 5 !< total number of gaseous tracers:
122                                            !< 1 = H2SO4, 2 = HNO3, 3 = NH3,
123                                            !< 4 = OCNV (non-volatile OC),
124                                            !< 5 = OCSV (semi-volatile) 
125    INTEGER(iwp), PARAMETER ::  nmod    = 7 !< number of modes for initialising
126                                            !< the aerosol size distribution                                             
127    INTEGER(iwp), PARAMETER ::  nreg    = 2 !< Number of main size subranges
128    INTEGER(iwp), PARAMETER ::  maxspec = 7 !< Max. number of aerosol species
129!   
130!-- Universal constants
131    REAL(wp), PARAMETER ::  abo    = 1.380662E-23_wp  !< Boltzmann constant (J/K)
132    REAL(wp), PARAMETER ::  alv    = 2.260E+6_wp      !< latent heat for H2O
133                                                      !< vaporisation (J/kg)
134    REAL(wp), PARAMETER ::  alv_d_rv  = 4896.96865_wp !< alv / rv
135    REAL(wp), PARAMETER ::  am_airmol = 4.8096E-26_wp !< Average mass of one air
136                                                      !< molecule (Jacobson,
137                                                      !< 2005, Eq. 2.3)                                                   
138    REAL(wp), PARAMETER ::  api6   = 0.5235988_wp     !< pi / 6   
139    REAL(wp), PARAMETER ::  argas  = 8.314409_wp      !< Gas constant (J/(mol K))
140    REAL(wp), PARAMETER ::  argas_d_cpd = 8.281283865E-3_wp !< argas per cpd
141    REAL(wp), PARAMETER ::  avo    = 6.02214E+23_wp   !< Avogadro constant (1/mol)
142    REAL(wp), PARAMETER ::  d_sa   = 5.539376964394570E-10_wp !< diameter of
143                                                      !< condensing sulphuric
144                                                      !< acid molecule (m) 
145    REAL(wp), PARAMETER ::  for_ppm_to_nconc =  7.243016311E+16_wp !<
146                                                 !< ppm * avo / R (K/(Pa*m3))
147    REAL(wp), PARAMETER ::  epsoc  = 0.15_wp          !< water uptake of organic
148                                                      !< material     
149    REAL(wp), PARAMETER ::  mclim  = 1.0E-23_wp    !< mass concentration min
150                                                   !< limit for aerosols (kg/m3)                                                   
151    REAL(wp), PARAMETER ::  n3     = 158.79_wp !< Number of H2SO4 molecules in
152                                               !< 3 nm cluster if d_sa=5.54e-10m
153    REAL(wp), PARAMETER ::  nclim  = 1.0_wp    !< number concentration min limit
154                                               !< for aerosols and gases (#/m3)
155    REAL(wp), PARAMETER ::  surfw0 = 0.073_wp  !< surface tension of pure water
156                                               !< at ~ 293 K (J/m2)   
157    REAL(wp), PARAMETER ::  vclim  = 1.0E-24_wp    !< volume concentration min
158                                                   !< limit for aerosols (m3/m3)                                           
159!-- Molar masses in kg/mol
160    REAL(wp), PARAMETER ::  ambc   = 12.0E-3_wp     !< black carbon (BC)
161    REAL(wp), PARAMETER ::  amdair = 28.970E-3_wp   !< dry air
162    REAL(wp), PARAMETER ::  amdu   = 100.E-3_wp     !< mineral dust
163    REAL(wp), PARAMETER ::  amh2o  = 18.0154E-3_wp  !< H2O
164    REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp  !< H2SO4
165    REAL(wp), PARAMETER ::  amhno3 = 63.01E-3_wp    !< HNO3
166    REAL(wp), PARAMETER ::  amn2o  = 44.013E-3_wp   !< N2O
167    REAL(wp), PARAMETER ::  amnh3  = 17.031E-3_wp   !< NH3
168    REAL(wp), PARAMETER ::  amo2   = 31.9988E-3_wp  !< O2
169    REAL(wp), PARAMETER ::  amo3   = 47.998E-3_wp   !< O3
170    REAL(wp), PARAMETER ::  amoc   = 150.E-3_wp     !< organic carbon (OC)
171    REAL(wp), PARAMETER ::  amss   = 58.44E-3_wp    !< sea salt (NaCl)
172!-- Densities in kg/m3
173    REAL(wp), PARAMETER ::  arhobc     = 2000.0_wp !< black carbon
174    REAL(wp), PARAMETER ::  arhodu     = 2650.0_wp !< mineral dust
175    REAL(wp), PARAMETER ::  arhoh2o    = 1000.0_wp !< H2O
176    REAL(wp), PARAMETER ::  arhoh2so4  = 1830.0_wp !< SO4
177    REAL(wp), PARAMETER ::  arhohno3   = 1479.0_wp !< HNO3
178    REAL(wp), PARAMETER ::  arhonh3    = 1530.0_wp !< NH3
179    REAL(wp), PARAMETER ::  arhooc     = 2000.0_wp !< organic carbon
180    REAL(wp), PARAMETER ::  arhoss     = 2165.0_wp !< sea salt (NaCl)
181!-- Volume of molecule in m3/#
182    REAL(wp), PARAMETER ::  amvh2o   = amh2o /avo / arhoh2o      !< H2O
183    REAL(wp), PARAMETER ::  amvh2so4 = amh2so4 / avo / arhoh2so4 !< SO4
184    REAL(wp), PARAMETER ::  amvhno3  = amhno3 / avo / arhohno3   !< HNO3
185    REAL(wp), PARAMETER ::  amvnh3   = amnh3 / avo / arhonh3     !< NH3 
186    REAL(wp), PARAMETER ::  amvoc    = amoc / avo / arhooc       !< OC
187    REAL(wp), PARAMETER ::  amvss    = amss / avo / arhoss       !< sea salt
188   
189!
190!-- SALSA switches:
191    INTEGER(iwp) ::  nj3 = 1 !< J3 parametrization (nucleation)
192                             !< 1 = condensational sink (Kerminen&Kulmala, 2002)
193                             !< 2 = coagulational sink (Lehtinen et al. 2007)
194                             !< 3 = coagS+self-coagulation (Anttila et al. 2010)                                       
195    INTEGER(iwp) ::  nsnucl = 0 !< Choice of the nucleation scheme:
196                                !< 0 = off   
197                                !< 1 = binary nucleation
198                                !< 2 = activation type nucleation
199                                !< 3 = kinetic nucleation
200                                !< 4 = ternary nucleation
201                                !< 5 = nucleation with ORGANICs
202                                !< 6 = activation type of nucleation with
203                                !<     H2SO4+ORG
204                                !< 7 = heteromolecular nucleation with H2SO4*ORG
205                                !< 8 = homomolecular nucleation of  H2SO4 +
206                                !<     heteromolecular nucleation with H2SO4*ORG
207                                !< 9 = homomolecular nucleation of  H2SO4 and ORG
208                                !<     +heteromolecular nucleation with H2SO4*ORG
209    LOGICAL ::  advect_particle_water = .TRUE.  !< advect water concentration of
210                                                !< particles                               
211    LOGICAL ::  decycle_lr            = .FALSE. !< Undo cyclic boundary
212                                                !< conditions: left and right
213    LOGICAL ::  decycle_ns            = .FALSE. !< north and south boundaries
214    LOGICAL ::  feedback_to_palm      = .FALSE. !< allow feedback due to
215                                                !< hydration and/or condensation
216                                                !< of H20
217    LOGICAL ::  no_insoluble          = .FALSE. !< Switch to exclude insoluble 
218                                                !< chemical components
219    LOGICAL ::  read_restart_data_salsa = .FALSE. !< read restart data for salsa
220    LOGICAL ::  salsa_gases_from_chem = .FALSE.   !< Transfer the gaseous
221                                                  !< components to SALSA from 
222                                                  !< from chemistry model
223    LOGICAL ::  van_der_waals_coagc   = .FALSE.   !< Enhancement of coagulation
224                                                  !< kernel by van der Waals and
225                                                  !< viscous forces
226    LOGICAL ::  write_binary_salsa    = .FALSE.   !< read binary for salsa
227!-- Process switches: nl* is read from the NAMELIST and is NOT changed.
228!--                   ls* is the switch used and will get the value of nl*
229!--                       except for special circumstances (spinup period etc.)
230    LOGICAL ::  nlcoag       = .FALSE. !< Coagulation master switch
231    LOGICAL ::  lscoag       = .FALSE. !<
232    LOGICAL ::  nlcnd        = .FALSE. !< Condensation master switch
233    LOGICAL ::  lscnd        = .FALSE. !<
234    LOGICAL ::  nlcndgas     = .FALSE. !< Condensation of precursor gases
235    LOGICAL ::  lscndgas     = .FALSE. !<
236    LOGICAL ::  nlcndh2oae   = .FALSE. !< Condensation of H2O on aerosol
237    LOGICAL ::  lscndh2oae   = .FALSE. !< particles (FALSE -> equilibrium calc.)
238    LOGICAL ::  nldepo       = .FALSE. !< Deposition master switch
239    LOGICAL ::  lsdepo       = .FALSE. !<
240    LOGICAL ::  nldepo_topo  = .FALSE. !< Deposition on vegetation master switch
241    LOGICAL ::  lsdepo_topo  = .FALSE. !<
242    LOGICAL ::  nldepo_vege  = .FALSE. !< Deposition on walls master switch
243    LOGICAL ::  lsdepo_vege  = .FALSE. !<
244    LOGICAL ::  nldistupdate = .TRUE.  !< Size distribution update master switch                                     
245    LOGICAL ::  lsdistupdate = .FALSE. !<                                     
246!
247!-- SALSA variables:
248    CHARACTER (LEN=20) ::  bc_salsa_b = 'neumann'   !< bottom boundary condition                                     
249    CHARACTER (LEN=20) ::  bc_salsa_t = 'neumann'   !< top boundary condition
250    CHARACTER (LEN=20) ::  depo_vege_type = 'zhang2001' !< or 'petroff2010'
251    CHARACTER (LEN=20) ::  depo_topo_type = 'zhang2001' !< or 'petroff2010'
252    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = & 
253                             (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
254                                 !< Decycling method at horizontal boundaries,
255                                 !< 1=left, 2=right, 3=south, 4=north
256                                 !< dirichlet = initial size distribution and
257                                 !< chemical composition set for the ghost and
258                                 !< first three layers
259                                 !< neumann = zero gradient
260    CHARACTER (LEN=3), DIMENSION(maxspec) ::  listspec = &  !< Active aerosols
261                                   (/'SO4','   ','   ','   ','   ','   ','   '/)
262    CHARACTER (LEN=20) ::  salsa_source_mode = 'no_source' 
263                                                    !< 'read_from_file',
264                                                    !< 'constant' or 'no_source'                                   
265    INTEGER(iwp) ::  dots_salsa = 0  !< starting index for salsa-timeseries
266    INTEGER(iwp) ::  fn1a = 1    !< last index for bin subranges:  subrange 1a
267    INTEGER(iwp) ::  fn2a = 1    !<                              subrange 2a
268    INTEGER(iwp) ::  fn2b = 1    !<                              subrange 2b
269    INTEGER(iwp), DIMENSION(ngast) ::  gas_index_chem = (/ 1, 1, 1, 1, 1/) !<
270                                 !< Index of gaseous compounds in the chemistry
271                                 !< model. In SALSA, 1 = H2SO4, 2 = HNO3,
272                                 !< 3 = NH3, 4 = OCNV, 5 = OCSV
273    INTEGER(iwp) ::  ibc_salsa_b !<
274    INTEGER(iwp) ::  ibc_salsa_t !<
275    INTEGER(iwp) ::  igctyp = 0  !< Initial gas concentration type
276                                 !< 0 = uniform (use H2SO4_init, HNO3_init,
277                                 !<     NH3_init, OCNV_init and OCSV_init)
278                                 !< 1 = read vertical profile from an input file 
279    INTEGER(iwp) ::  in1a = 1    !< start index for bin subranges: subrange 1a
280    INTEGER(iwp) ::  in2a = 1    !<                              subrange 2a
281    INTEGER(iwp) ::  in2b = 1    !<                              subrange 2b
282    INTEGER(iwp) ::  isdtyp = 0  !< Initial size distribution type
283                                 !< 0 = uniform
284                                 !< 1 = read vertical profile of the mode number
285                                 !<     concentration from an input file 
286    INTEGER(iwp) ::  ibc  = -1 !< Indice for: black carbon (BC)
287    INTEGER(iwp) ::  idu  = -1 !< dust
288    INTEGER(iwp) ::  inh  = -1 !< NH3
289    INTEGER(iwp) ::  ino  = -1 !< HNO3   
290    INTEGER(iwp) ::  ioc  = -1 !< organic carbon (OC)
291    INTEGER(iwp) ::  iso4 = -1 !< SO4 or H2SO4   
292    INTEGER(iwp) ::  iss  = -1 !< sea salt
293    INTEGER(iwp) ::  lod_aero = 0   !< level of detail for aerosol emissions
294    INTEGER(iwp) ::  lod_gases = 0  !< level of detail for gaseous emissions   
295    INTEGER(iwp), DIMENSION(nreg) ::  nbin = (/ 3, 7/)    !< Number of size bins
296                                               !< for each aerosol size subrange
297    INTEGER(iwp) ::  nbins = 1  !< total number of size bins
298    INTEGER(iwp) ::  ncc   = 1  !< number of chemical components used     
299    INTEGER(iwp) ::  ncc_tot = 1!< total number of chemical compounds (ncc+1
300                                !< if particle water is advected)
301    REAL(wp) ::  act_coeff = 1.0E-7_wp     !< Activation coefficient
302    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  emission_mass_fracs  !< array for
303                                    !< aerosol composition per emission category
304                                    !< 1:SO4 2:OC 3:BC 4:DU 5:SS 6:NO 7:NH 
305    REAL(wp) ::  dt_salsa  = 0.00001_wp    !< Time step of SALSA
306    REAL(wp) ::  H2SO4_init = nclim        !< Init value for sulphuric acid gas
307    REAL(wp) ::  HNO3_init  = nclim        !< Init value for nitric acid gas
308    REAL(wp) ::  last_salsa_time = 0.0_wp  !< time of the previous salsa
309                                           !< timestep
310    REAL(wp) ::  nf2a = 1.0_wp             !< Number fraction allocated to a-
311                                           !< bins in subrange 2
312                                           !< (b-bins will get 1-nf2a)   
313    REAL(wp) ::  NH3_init  = nclim         !< Init value for ammonia gas
314    REAL(wp) ::  OCNV_init = nclim         !< Init value for non-volatile
315                                           !< organic gases
316    REAL(wp) ::  OCSV_init = nclim         !< Init value for semi-volatile
317                                           !< organic gases
318    REAL(wp), DIMENSION(nreg+1) ::  reglim = & !< Min&max diameters of size subranges
319                                 (/ 3.0E-9_wp, 5.0E-8_wp, 1.0E-5_wp/)
320    REAL(wp) ::  rhlim = 1.20_wp    !< RH limit in %/100. Prevents
321                                    !< unrealistically high RH in condensation                           
322    REAL(wp) ::  skip_time_do_salsa = 0.0_wp !< Starting time of SALSA (s)
323!-- Initial log-normal size distribution: mode diameter (dpg, micrometres),
324!-- standard deviation (sigmag) and concentration (n_lognorm, #/cm3)
325    REAL(wp), DIMENSION(nmod) ::  dpg   = (/0.013_wp, 0.054_wp, 0.86_wp,       &
326                                            0.2_wp, 0.2_wp, 0.2_wp, 0.2_wp/) 
327    REAL(wp), DIMENSION(nmod) ::  sigmag  = (/1.8_wp, 2.16_wp, 2.21_wp,        &
328                                              2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/) 
329    REAL(wp), DIMENSION(nmod) ::  n_lognorm = (/1.04e+5_wp, 3.23E+4_wp, 5.4_wp,&
330                                                0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
331!-- Initial mass fractions / chemical composition of the size distribution   
332    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_a = & !< mass fractions between
333             (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for A bins
334    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_b = & !< mass fractions between
335             (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for B bins
336             
337    REAL(wp), ALLOCATABLE, DIMENSION(:) ::  bin_low_limits  !< to deliver
338                                                            !< information about
339                                                            !< the lower
340                                                            !< diameters per bin                                       
341    REAL(wp), ALLOCATABLE, DIMENSION(:) ::  nsect     !< Background number
342                                                      !< concentration per bin
343    REAL(wp), ALLOCATABLE, DIMENSION(:) ::  massacc   !< Mass accomodation
344                                                      !< coefficients per bin                                             
345!
346!-- SALSA derived datatypes:
347!
348!-- Prognostic variable: Aerosol size bin information (number (#/m3) and
349!-- mass (kg/m3) concentration) and the concentration of gaseous tracers (#/m3).
350!-- Gas tracers are contained sequentially in dimension 4 as:
351!-- 1. H2SO4, 2. HNO3, 3. NH3, 4. OCNV (non-volatile organics),
352!-- 5. OCSV (semi-volatile)
353    TYPE salsa_variable
354       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS     ::  conc
355       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS     ::  conc_p
356       REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS     ::  tconc_m
357       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  flux_s, diss_s
358       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l, diss_l
359       REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  init
360       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  source
361       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  sums_ws_l
362    END TYPE salsa_variable
363   
364!-- Map bin indices between parallel size distributions   
365    TYPE t_parallelbin
366       INTEGER(iwp) ::  cur  ! Index for current distribution
367       INTEGER(iwp) ::  par  ! Index for corresponding parallel distribution
368    END TYPE t_parallelbin
369   
370!-- Datatype used to store information about the binned size distributions of
371!-- aerosols
372    TYPE t_section
373       REAL(wp) ::  vhilim   !< bin volume at the high limit
374       REAL(wp) ::  vlolim   !< bin volume at the low limit
375       REAL(wp) ::  vratiohi !< volume ratio between the center and high limit
376       REAL(wp) ::  vratiolo !< volume ratio between the center and low limit
377       REAL(wp) ::  dmid     !< bin middle diameter (m)
378       !******************************************************
379       ! ^ Do NOT change the stuff above after initialization !
380       !******************************************************
381       REAL(wp) ::  dwet    !< Wet diameter or mean droplet diameter (m)
382       REAL(wp), DIMENSION(maxspec+1) ::  volc !< Volume concentrations
383                            !< (m^3/m^3) of aerosols + water. Since most of
384                            !< the stuff in SALSA is hard coded, these *have to
385                            !< be* in the order
386                            !< 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
387       REAL(wp) ::  veqh2o  !< Equilibrium H2O concentration for each particle
388       REAL(wp) ::  numc    !< Number concentration of particles/droplets (#/m3)
389       REAL(wp) ::  core    !< Volume of dry particle
390    END TYPE t_section 
391!
392!-- Local aerosol properties in SALSA
393    TYPE(t_section), ALLOCATABLE ::  aero(:)
394!
395!-- SALSA tracers:
396!-- Tracers as x = x(k,j,i,bin). The 4th dimension contains all the size bins
397!-- sequentially for each aerosol species  + water.
398!
399!-- Prognostic tracers:
400!
401!-- Number concentration (#/m3)
402    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_number
403    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_1
404    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_2
405    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_3
406!
407!-- Mass concentration (kg/m3)
408    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_mass
409    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_1
410    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_2
411    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_3
412!
413!-- Gaseous tracers (#/m3)
414    TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  salsa_gas
415    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_1
416    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_2
417    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_3
418!
419!-- Diagnostic tracers
420    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  sedim_vd !< sedimentation
421                                                           !< velocity per size
422                                                           !< bin (m/s)
423    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  Ra_dry !< dry radius (m)
424   
425!-- Particle component index tables
426    TYPE(component_index) :: prtcl !< Contains "getIndex" which gives the index
427                                   !< for a given aerosol component name, i.e.
428                                   !< 1:SO4, 2:OC, 3:BC, 4:DU,
429                                   !< 5:SS, 6:NO, 7:NH, 8:H2O 
430!                                   
431!-- Data output arrays:
432!-- Gases:
433    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_H2SO4_av  !< H2SO4
434    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_HNO3_av   !< HNO3
435    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_NH3_av    !< NH3
436    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_OCNV_av   !< non-vola-
437                                                                    !< tile OC
438    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_OCSV_av   !< semi-vol.
439                                                                    !< OC
440!-- Integrated:                                                                   
441    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  LDSA_av  !< lung-
442                                                                 !< deposited
443                                                                 !< surface area                                                   
444    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  Ntot_av  !< total number
445                                                                 !< conc.
446    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  PM25_av  !< PM2.5
447    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  PM10_av  !< PM10
448!-- In the particle phase:   
449    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_BC_av  !< black carbon
450    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_DU_av  !< dust
451    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_H2O_av !< liquid water
452    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_NH_av  !< ammonia
453    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_NO_av  !< nitrates
454    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_OC_av  !< org. carbon
455    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_SO4_av !< sulphates
456    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_SS_av  !< sea salt
457!-- Bins:   
458    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mbins_av !< bin mass 
459    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  Nbins_av !< bin number
460
461   
462!
463!-- PALM interfaces:
464!
465!-- Boundary conditions:
466    INTERFACE salsa_boundary_conds
467       MODULE PROCEDURE salsa_boundary_conds
468       MODULE PROCEDURE salsa_boundary_conds_decycle
469    END INTERFACE salsa_boundary_conds
470!   
471!-- Data output checks for 2D/3D data to be done in check_parameters
472    INTERFACE salsa_check_data_output
473       MODULE PROCEDURE salsa_check_data_output
474    END INTERFACE salsa_check_data_output
475   
476!
477!-- Input parameter checks to be done in check_parameters
478    INTERFACE salsa_check_parameters
479       MODULE PROCEDURE salsa_check_parameters
480    END INTERFACE salsa_check_parameters
481
482!
483!-- Averaging of 3D data for output
484    INTERFACE salsa_3d_data_averaging
485       MODULE PROCEDURE salsa_3d_data_averaging
486    END INTERFACE salsa_3d_data_averaging
487
488!
489!-- Data output of 2D quantities
490    INTERFACE salsa_data_output_2d
491       MODULE PROCEDURE salsa_data_output_2d
492    END INTERFACE salsa_data_output_2d
493
494!
495!-- Data output of 3D data
496    INTERFACE salsa_data_output_3d
497       MODULE PROCEDURE salsa_data_output_3d
498    END INTERFACE salsa_data_output_3d
499   
500!
501!-- Data output of 3D data
502    INTERFACE salsa_data_output_mask
503       MODULE PROCEDURE salsa_data_output_mask
504    END INTERFACE salsa_data_output_mask
505
506!
507!-- Definition of data output quantities
508    INTERFACE salsa_define_netcdf_grid
509       MODULE PROCEDURE salsa_define_netcdf_grid
510    END INTERFACE salsa_define_netcdf_grid
511   
512!
513!-- Output of information to the header file
514    INTERFACE salsa_header
515       MODULE PROCEDURE salsa_header
516    END INTERFACE salsa_header
517 
518!
519!-- Initialization actions 
520    INTERFACE salsa_init
521       MODULE PROCEDURE salsa_init
522    END INTERFACE salsa_init
523 
524!
525!-- Initialization of arrays
526    INTERFACE salsa_init_arrays
527       MODULE PROCEDURE salsa_init_arrays
528    END INTERFACE salsa_init_arrays
529
530!
531!-- Writing of binary output for restart runs  !!! renaming?!
532    INTERFACE salsa_wrd_local
533       MODULE PROCEDURE salsa_wrd_local
534    END INTERFACE salsa_wrd_local
535   
536!
537!-- Reading of NAMELIST parameters
538    INTERFACE salsa_parin
539       MODULE PROCEDURE salsa_parin
540    END INTERFACE salsa_parin
541
542!
543!-- Reading of parameters for restart runs
544    INTERFACE salsa_rrd_local
545       MODULE PROCEDURE salsa_rrd_local
546    END INTERFACE salsa_rrd_local
547   
548!
549!-- Swapping of time levels (required for prognostic variables)
550    INTERFACE salsa_swap_timelevel
551       MODULE PROCEDURE salsa_swap_timelevel
552    END INTERFACE salsa_swap_timelevel
553
554    INTERFACE salsa_driver
555       MODULE PROCEDURE salsa_driver
556    END INTERFACE salsa_driver
557
558    INTERFACE salsa_tendency
559       MODULE PROCEDURE salsa_tendency
560       MODULE PROCEDURE salsa_tendency_ij
561    END INTERFACE salsa_tendency
562   
563   
564   
565    SAVE
566
567    PRIVATE
568!
569!-- Public functions:
570    PUBLIC salsa_boundary_conds, salsa_check_data_output,                      &
571           salsa_check_parameters, salsa_3d_data_averaging,                    &
572           salsa_data_output_2d, salsa_data_output_3d, salsa_data_output_mask, &
573           salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver,          &
574           salsa_header, salsa_init, salsa_init_arrays, salsa_parin,           &
575           salsa_rrd_local, salsa_swap_timelevel, salsa_tendency,              &
576           salsa_wrd_local
577!
578!-- Public parameters, constants and initial values
579    PUBLIC dots_salsa, dt_salsa, last_salsa_time, lsdepo, salsa,               &
580           salsa_gases_from_chem, skip_time_do_salsa
581!
582!-- Public prognostic variables
583    PUBLIC aerosol_mass, aerosol_number, fn2a, fn2b, gconc_2, in1a, in2b,      &
584           mconc_2, nbins, ncc, ncc_tot, nclim, nconc_2, ngast, prtcl, Ra_dry, &
585           salsa_gas, sedim_vd
586           
587
588 CONTAINS
589
590!------------------------------------------------------------------------------!
591! Description:
592! ------------
593!> Parin for &salsa_par for new modules
594!------------------------------------------------------------------------------!
595 SUBROUTINE salsa_parin
596
597    IMPLICIT NONE
598
599    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line
600                                  !< of the parameter file
601                                 
602    NAMELIST /salsa_parameters/             &
603                          advect_particle_water, & ! Switch for advecting
604                                                ! particle water. If .FALSE.,
605                                                ! equilibration is called at
606                                                ! each time step.       
607                          bc_salsa_b,       &   ! bottom boundary condition
608                          bc_salsa_t,       &   ! top boundary condition
609                          decycle_lr,       &   ! decycle SALSA components
610                          decycle_method,   &   ! decycle method applied:
611                                                ! 1=left 2=right 3=south 4=north
612                          decycle_ns,       &   ! decycle SALSA components
613                          depo_vege_type,   &   ! Parametrisation type
614                          depo_topo_type,   &   ! Parametrisation type
615                          dpg,              &   ! Mean diameter for the initial
616                                                ! log-normal modes
617                          dt_salsa,         &   ! SALSA timestep in seconds
618                          feedback_to_palm, &   ! allow feedback due to
619                                                ! hydration / condensation
620                          H2SO4_init,       &   ! Init value for sulphuric acid
621                          HNO3_init,        &   ! Init value for nitric acid
622                          igctyp,           &   ! Initial gas concentration type
623                          isdtyp,           &   ! Initial size distribution type                                               
624                          listspec,         &   ! List of actived aerosols
625                                                ! (string list)
626                          mass_fracs_a,     &   ! Initial relative contribution 
627                                                ! of each species to particle 
628                                                ! volume in a-bins, 0 for unused
629                          mass_fracs_b,     &   ! Initial relative contribution 
630                                                ! of each species to particle
631                                                ! volume in b-bins, 0 for unused
632                          n_lognorm,        &   ! Number concentration for the
633                                                ! log-normal modes                                               
634                          nbin,             &   ! Number of size bins for
635                                                ! aerosol size subranges 1 & 2
636                          nf2a,             &   ! Number fraction of particles
637                                                ! allocated to a-bins in
638                                                ! subrange 2 b-bins will get
639                                                ! 1-nf2a                         
640                          NH3_init,         &   ! Init value for ammonia
641                          nj3,              &   ! J3 parametrization
642                                                ! 1 = condensational sink
643                                                !     (Kerminen&Kulmala, 2002)
644                                                ! 2 = coagulational sink
645                                                !     (Lehtinen et al. 2007)
646                                                ! 3 = coagS+self-coagulation
647                                                !     (Anttila et al. 2010)                                                   
648                          nlcnd,            &   ! Condensation master switch
649                          nlcndgas,         &   ! Condensation of gases
650                          nlcndh2oae,       &   ! Condensation of H2O                           
651                          nlcoag,           &   ! Coagulation master switch
652                          nldepo,           &   ! Deposition master switch
653                          nldepo_vege,      &   ! Deposition on vegetation
654                                                ! master switch
655                          nldepo_topo,      &   ! Deposition on topo master
656                                                ! switch                         
657                          nldistupdate,     &   ! Size distribution update
658                                                ! master switch
659                          nsnucl,           &   ! Nucleation scheme:
660                                                ! 0 = off,
661                                                ! 1 = binary nucleation
662                                                ! 2 = activation type nucleation
663                                                ! 3 = kinetic nucleation
664                                                ! 4 = ternary nucleation
665                                                ! 5 = nucleation with organics
666                                                ! 6 = activation type of
667                                                !     nucleation with H2SO4+ORG
668                                                ! 7 = heteromolecular nucleation
669                                                !     with H2SO4*ORG
670                                                ! 8 = homomolecular nucleation 
671                                                !     of H2SO4 + heteromolecular
672                                                !     nucleation with H2SO4*ORG
673                                                ! 9 = homomolecular nucleation
674                                                !     of H2SO4 and ORG + hetero-
675                                                !     molecular nucleation with
676                                                !     H2SO4*ORG
677                          OCNV_init,        &   ! Init value for non-volatile
678                                                ! organic gases
679                          OCSV_init,        &   ! Init value for semi-volatile
680                                                ! organic gases
681                          read_restart_data_salsa, & ! read restart data for
682                                                     ! salsa
683                          reglim,           &   ! Min&max diameter limits of
684                                                ! size subranges
685                          salsa,            &   ! Master switch for SALSA
686                          salsa_source_mode,&   ! 'read_from_file' or 'constant'
687                                                ! or 'no_source'
688                          sigmag,           &   ! stdev for the initial log-
689                                                ! normal modes                                               
690                          skip_time_do_salsa, & ! Starting time of SALSA (s)
691                          van_der_waals_coagc,& ! include van der Waals forces
692                          write_binary_salsa    ! Write binary for salsa
693                           
694       
695    line = ' '
696       
697!
698!-- Try to find salsa package
699    REWIND ( 11 )
700    line = ' '
701    DO WHILE ( INDEX( line, '&salsa_parameters' ) == 0 )
702       READ ( 11, '(A)', END=10 )  line
703    ENDDO
704    BACKSPACE ( 11 )
705
706!
707!-- Read user-defined namelist
708    READ ( 11, salsa_parameters )
709
710!
711!-- Enable salsa (salsa switch in modules.f90)
712    salsa = .TRUE.
713
714 10 CONTINUE
715       
716 END SUBROUTINE salsa_parin
717
718 
719!------------------------------------------------------------------------------!
720! Description:
721! ------------
722!> Check parameters routine for salsa.
723!------------------------------------------------------------------------------!
724 SUBROUTINE salsa_check_parameters
725
726    USE control_parameters,                                                    &
727        ONLY:  message_string
728       
729    IMPLICIT NONE
730   
731!
732!-- Checks go here (cf. check_parameters.f90).
733    IF ( salsa  .AND.  .NOT.  humidity )  THEN
734       WRITE( message_string, * ) 'salsa = ', salsa, ' is ',                   &
735              'not allowed with humidity = ', humidity
736       CALL message( 'check_parameters', 'SA0009', 1, 2, 0, 6, 0 )
737    ENDIF
738   
739    IF ( bc_salsa_b == 'dirichlet' )  THEN
740       ibc_salsa_b = 0
741    ELSEIF ( bc_salsa_b == 'neumann' )  THEN
742       ibc_salsa_b = 1
743    ELSE
744       message_string = 'unknown boundary condition: bc_salsa_b = "'           &
745                         // TRIM( bc_salsa_t ) // '"'
746       CALL message( 'check_parameters', 'SA0011', 1, 2, 0, 6, 0 )                 
747    ENDIF
748   
749    IF ( bc_salsa_t == 'dirichlet' )  THEN
750       ibc_salsa_t = 0
751    ELSEIF ( bc_salsa_t == 'neumann' )  THEN
752       ibc_salsa_t = 1
753    ELSE
754       message_string = 'unknown boundary condition: bc_salsa_t = "'           &
755                         // TRIM( bc_salsa_t ) // '"'
756       CALL message( 'check_parameters', 'SA0012', 1, 2, 0, 6, 0 )                 
757    ENDIF
758   
759    IF ( nj3 < 1  .OR.  nj3 > 3 )  THEN
760       message_string = 'unknown nj3 (must be 1-3)'
761       CALL message( 'check_parameters', 'SA0044', 1, 2, 0, 6, 0 )
762    ENDIF
763           
764 END SUBROUTINE salsa_check_parameters
765
766!------------------------------------------------------------------------------!
767!
768! Description:
769! ------------
770!> Subroutine defining appropriate grid for netcdf variables.
771!> It is called out from subroutine netcdf.
772!> Same grid as for other scalars (see netcdf_interface_mod.f90)
773!------------------------------------------------------------------------------!
774 SUBROUTINE salsa_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
775   
776    IMPLICIT NONE
777
778    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !<
779    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !<
780    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !<
781    CHARACTER (LEN=*), INTENT(IN)  ::  var      !<
782   
783    LOGICAL, INTENT(OUT) ::  found   !<
784   
785    found  = .TRUE.
786!
787!-- Check for the grid
788
789    IF ( var(1:2) == 'g_' )  THEN
790       grid_x = 'x' 
791       grid_y = 'y' 
792       grid_z = 'zu'   
793    ELSEIF ( var(1:4) == 'LDSA' )  THEN
794       grid_x = 'x' 
795       grid_y = 'y' 
796       grid_z = 'zu'
797    ELSEIF ( var(1:5) == 'm_bin' )  THEN
798       grid_x = 'x' 
799       grid_y = 'y' 
800       grid_z = 'zu'
801    ELSEIF ( var(1:5) == 'N_bin' )  THEN
802       grid_x = 'x' 
803       grid_y = 'y' 
804       grid_z = 'zu'
805    ELSEIF ( var(1:4) == 'Ntot' ) THEN
806       grid_x = 'x' 
807       grid_y = 'y' 
808       grid_z = 'zu'
809    ELSEIF ( var(1:2) == 'PM' )  THEN
810       grid_x = 'x' 
811       grid_y = 'y' 
812       grid_z = 'zu'
813    ELSEIF ( var(1:2) == 's_' )  THEN
814       grid_x = 'x' 
815       grid_y = 'y' 
816       grid_z = 'zu'
817    ELSE
818       found  = .FALSE.
819       grid_x = 'none'
820       grid_y = 'none'
821       grid_z = 'none'
822    ENDIF
823
824 END SUBROUTINE salsa_define_netcdf_grid
825
826 
827!------------------------------------------------------------------------------!
828! Description:
829! ------------
830!> Header output for new module
831!------------------------------------------------------------------------------!
832 SUBROUTINE salsa_header( io )
833
834    IMPLICIT NONE
835 
836    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
837!
838!-- Write SALSA header
839    WRITE( io, 1 )
840    WRITE( io, 2 ) skip_time_do_salsa
841    WRITE( io, 3 ) dt_salsa
842    WRITE( io, 12 )  SHAPE( aerosol_number(1)%conc ), nbins
843    IF ( advect_particle_water )  THEN
844       WRITE( io, 16 )  SHAPE( aerosol_mass(1)%conc ), ncc_tot*nbins,          &
845                        advect_particle_water
846    ELSE
847       WRITE( io, 16 )  SHAPE( aerosol_mass(1)%conc ), ncc*nbins,              &
848                        advect_particle_water
849    ENDIF
850    IF ( .NOT. salsa_gases_from_chem )  THEN
851       WRITE( io, 17 )  SHAPE( aerosol_mass(1)%conc ), ngast,                  &
852                        salsa_gases_from_chem
853    ENDIF
854    WRITE( io, 4 ) 
855    IF ( nsnucl > 0 )  THEN
856       WRITE( io, 5 ) nsnucl, nj3
857    ENDIF
858    IF ( nlcoag )  THEN
859       WRITE( io, 6 ) 
860    ENDIF
861    IF ( nlcnd )  THEN
862       WRITE( io, 7 ) nlcndgas, nlcndh2oae
863    ENDIF
864    IF ( nldepo )  THEN
865       WRITE( io, 14 ) nldepo_vege, nldepo_topo
866    ENDIF
867    WRITE( io, 8 )  reglim, nbin, bin_low_limits
868    WRITE( io, 15 ) nsect
869    WRITE( io, 13 ) ncc, listspec, mass_fracs_a, mass_fracs_b
870    IF ( .NOT. salsa_gases_from_chem )  THEN
871       WRITE( io, 18 ) ngast, H2SO4_init, HNO3_init, NH3_init, OCNV_init,      &
872                       OCSV_init
873    ENDIF
874    WRITE( io, 9 )  isdtyp, igctyp
875    IF ( isdtyp == 0 )  THEN
876       WRITE( io, 10 )  dpg, sigmag, n_lognorm
877    ELSE
878       WRITE( io, 11 )
879    ENDIF
880   
881
8821   FORMAT (//' SALSA information:'/                                           &
883              ' ------------------------------'/)
8842   FORMAT   ('    Starts at: skip_time_do_salsa = ', F10.2, '  s')
8853   FORMAT  (/'    Timestep: dt_salsa = ', F6.2, '  s')
88612  FORMAT  (/'    Array shape (z,y,x,bins):'/                                 &
887              '       aerosol_number:  ', 4(I3)) 
88816  FORMAT  (/'       aerosol_mass:    ', 4(I3),/                              &
889              '       (advect_particle_water = ', L1, ')')
89017  FORMAT   ('       salsa_gas: ', 4(I3),/                                    &
891              '       (salsa_gases_from_chem = ', L1, ')')
8924   FORMAT  (/'    Aerosol dynamic processes included: ')
8935   FORMAT  (/'       nucleation (scheme = ', I1, ' and J3 parametrization = ',&
894               I1, ')')
8956   FORMAT  (/'       coagulation')
8967   FORMAT  (/'       condensation (of precursor gases = ', L1,                &
897              '          and water vapour = ', L1, ')' )
89814  FORMAT  (/'       dry deposition (on vegetation = ', L1,                   &
899              '          and on topography = ', L1, ')')             
9008   FORMAT  (/'    Aerosol bin subrange limits (in metres): ',  3(ES10.2E3), / &
901              '    Number of size bins for each aerosol subrange: ', 2I3,/     &
902              '    Aerosol bin limits (in metres): ', 9(ES10.2E3))
90315  FORMAT   ('    Initial number concentration in bins at the lowest level',  &
904              ' (#/m**3):', 9(ES10.2E3))       
90513  FORMAT  (/'    Number of chemical components used: ', I1,/                 &
906              '       Species: ',7(A6),/                                       &
907              '    Initial relative contribution of each species to particle', & 
908              ' volume in:',/                                                  &
909              '       a-bins: ', 7(F6.3),/                                     &
910              '       b-bins: ', 7(F6.3))
91118  FORMAT  (/'    Number of gaseous tracers used: ', I1,/                     &
912              '    Initial gas concentrations:',/                              &
913              '       H2SO4: ',ES12.4E3, ' #/m**3',/                           &
914              '       HNO3:  ',ES12.4E3, ' #/m**3',/                           &
915              '       NH3:   ',ES12.4E3, ' #/m**3',/                           &
916              '       OCNV:  ',ES12.4E3, ' #/m**3',/                           &
917              '       OCSV:  ',ES12.4E3, ' #/m**3')
9189    FORMAT (/'   Initialising concentrations: ', /                            &
919              '      Aerosol size distribution: isdtyp = ', I1,/               &
920              '      Gas concentrations: igctyp = ', I1 )
92110   FORMAT ( '      Mode diametres: dpg(nmod) = ', 7(F7.3),/                  &
922              '      Standard deviation: sigmag(nmod) = ', 7(F7.2),/           &
923              '      Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3) )
92411   FORMAT (/'      Size distribution read from a file.')
925
926 END SUBROUTINE salsa_header
927
928!------------------------------------------------------------------------------!
929! Description:
930! ------------
931!> Allocate SALSA arrays and define pointers if required
932!------------------------------------------------------------------------------!
933 SUBROUTINE salsa_init_arrays
934 
935    USE surface_mod,                                                           &
936        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
937               surf_usm_v
938
939    IMPLICIT NONE
940   
941    INTEGER(iwp) ::  gases_available !< Number of available gas components in
942                                     !< the chemistry model
943    INTEGER(iwp) ::  i   !< loop index for allocating
944    INTEGER(iwp) ::  l   !< loop index for allocating: surfaces
945    INTEGER(iwp) ::  lsp !< loop index for chem species in the chemistry model
946   
947    gases_available = 0
948
949!
950!-- Allocate prognostic variables (see salsa_swap_timelevel)
951
952!
953!-- Set derived indices:
954!-- (This does the same as the subroutine salsa_initialize in SALSA/
955!-- UCLALES-SALSA)       
956    in1a = 1                ! 1st index of subrange 1a
957    in2a = in1a + nbin(1)   ! 1st index of subrange 2a
958    fn1a = in2a - 1         ! last index of subrange 1a
959    fn2a = fn1a + nbin(2)   ! last index of subrange 2a
960   
961!   
962!-- If the fraction of insoluble aerosols in subrange 2 is zero: do not allocate
963!-- arrays for them
964    IF ( nf2a > 0.999999_wp  .AND.  SUM( mass_fracs_b ) < 0.00001_wp )  THEN
965       no_insoluble = .TRUE.
966       in2b = fn2a+1    ! 1st index of subrange 2b
967       fn2b = fn2a      ! last index of subrange 2b
968    ELSE
969       in2b = in2a + nbin(2)   ! 1st index of subrange 2b
970       fn2b = fn2a + nbin(2)   ! last index of subrange 2b
971    ENDIF
972   
973   
974    nbins = fn2b   ! total number of aerosol size bins
975!   
976!-- Create index tables for different aerosol components
977    CALL component_index_constructor( prtcl, ncc, maxspec, listspec )
978   
979    ncc_tot = ncc
980    IF ( advect_particle_water )  ncc_tot = ncc + 1  ! Add water
981   
982!
983!-- Allocate:
984    ALLOCATE( aero(nbins), bin_low_limits(nbins), nsect(nbins), massacc(nbins) )
985    IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins) )         
986    ALLOCATE( Ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins) )
987   
988!   
989!-- Aerosol number concentration
990    ALLOCATE( aerosol_number(nbins) )
991    ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins),                    &
992              nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins),                    &
993              nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins) )
994    nconc_1 = 0.0_wp
995    nconc_2 = 0.0_wp
996    nconc_3 = 0.0_wp
997   
998    DO i = 1, nbins
999       aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => nconc_1(:,:,:,i)
1000       aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => nconc_2(:,:,:,i)
1001       aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i)
1002       ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),     &
1003                 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),     &
1004                 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1005                 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1006                 aerosol_number(i)%init(nzb:nzt+1),                            &
1007                 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1008    ENDDO     
1009   
1010!   
1011!-- Aerosol mass concentration   
1012    ALLOCATE( aerosol_mass(ncc_tot*nbins) ) 
1013    ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncc_tot*nbins),            &
1014              mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncc_tot*nbins),            &
1015              mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncc_tot*nbins) )
1016    mconc_1 = 0.0_wp
1017    mconc_2 = 0.0_wp
1018    mconc_3 = 0.0_wp
1019   
1020    DO i = 1, ncc_tot*nbins
1021       aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => mconc_1(:,:,:,i)
1022       aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => mconc_2(:,:,:,i)
1023       aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i)       
1024       ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
1025                 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
1026                 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1027                 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1028                 aerosol_mass(i)%init(nzb:nzt+1),                              &
1029                 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
1030    ENDDO
1031   
1032!
1033!-- Surface fluxes: answs = aerosol number, amsws = aerosol mass
1034!
1035!-- Horizontal surfaces: default type
1036    DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
1037       ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins ) )
1038       ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins*ncc_tot ) )
1039       surf_def_h(l)%answs = 0.0_wp
1040       surf_def_h(l)%amsws = 0.0_wp
1041    ENDDO
1042!-- Horizontal surfaces: natural type   
1043    IF ( land_surface )  THEN
1044       ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins ) )
1045       ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins*ncc_tot ) )
1046       surf_lsm_h%answs = 0.0_wp
1047       surf_lsm_h%amsws = 0.0_wp
1048    ENDIF
1049!-- Horizontal surfaces: urban type
1050    IF ( urban_surface )  THEN
1051       ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins ) )
1052       ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins*ncc_tot ) )
1053       surf_usm_h%answs = 0.0_wp
1054       surf_usm_h%amsws = 0.0_wp
1055    ENDIF
1056!
1057!-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
1058!-- westward (l=3) facing
1059    DO  l = 0, 3   
1060       ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins ) )
1061       surf_def_v(l)%answs = 0.0_wp
1062       ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins*ncc_tot ) )
1063       surf_def_v(l)%amsws = 0.0_wp
1064       
1065       IF ( land_surface)  THEN
1066          ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins ) )
1067          surf_lsm_v(l)%answs = 0.0_wp
1068          ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins*ncc_tot ) )
1069          surf_lsm_v(l)%amsws = 0.0_wp
1070       ENDIF
1071       
1072       IF ( urban_surface )  THEN
1073          ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins ) )
1074          surf_usm_v(l)%answs = 0.0_wp
1075          ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins*ncc_tot ) )
1076          surf_usm_v(l)%amsws = 0.0_wp
1077       ENDIF
1078    ENDDO   
1079   
1080!
1081!-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV)
1082!-- (number concentration (#/m3) )
1083!
1084!-- If chemistry is on, read gas phase concentrations from there. Otherwise,
1085!-- allocate salsa_gas array.
1086
1087    IF ( air_chemistry )  THEN   
1088       DO  lsp = 1, nvar
1089          IF ( TRIM( chem_species(lsp)%name ) == 'H2SO4' )  THEN
1090             gases_available = gases_available + 1
1091             gas_index_chem(1) = lsp
1092          ELSEIF ( TRIM( chem_species(lsp)%name ) == 'HNO3' )  THEN
1093             gases_available = gases_available + 1 
1094             gas_index_chem(2) = lsp
1095          ELSEIF ( TRIM( chem_species(lsp)%name ) == 'NH3' )  THEN
1096             gases_available = gases_available + 1
1097             gas_index_chem(3) = lsp
1098          ELSEIF ( TRIM( chem_species(lsp)%name ) == 'OCNV' )  THEN
1099             gases_available = gases_available + 1
1100             gas_index_chem(4) = lsp
1101          ELSEIF ( TRIM( chem_species(lsp)%name ) == 'OCSV' )  THEN
1102             gases_available = gases_available + 1
1103             gas_index_chem(5) = lsp
1104          ENDIF
1105       ENDDO
1106
1107       IF ( gases_available == ngast )  THEN
1108          salsa_gases_from_chem = .TRUE.
1109       ELSE
1110          WRITE( message_string, * ) 'SALSA is run together with chemistry '// &
1111                                     'but not all gaseous components are '//   &
1112                                     'provided by kpp (H2SO4, HNO3, NH3, '//   &
1113                                     'OCNV, OCSC)'
1114       CALL message( 'check_parameters', 'SA0024', 1, 2, 0, 6, 0 )
1115       ENDIF
1116
1117    ELSE
1118
1119       ALLOCATE( salsa_gas(ngast) ) 
1120       ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngast),                 &
1121                 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngast),                 &
1122                 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngast) )
1123       gconc_1 = 0.0_wp
1124       gconc_2 = 0.0_wp
1125       gconc_3 = 0.0_wp
1126       
1127       DO i = 1, ngast
1128          salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => gconc_1(:,:,:,i)
1129          salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => gconc_2(:,:,:,i)
1130          salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i)
1131          ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
1132                    salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
1133                    salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1134                    salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
1135                    salsa_gas(i)%init(nzb:nzt+1),                              &
1136                    salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1137       ENDDO       
1138!
1139!--    Surface fluxes: gtsws = gaseous tracer flux
1140!
1141!--    Horizontal surfaces: default type
1142       DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
1143          ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngast ) )
1144          surf_def_h(l)%gtsws = 0.0_wp
1145       ENDDO
1146!--    Horizontal surfaces: natural type   
1147       IF ( land_surface )  THEN
1148          ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngast ) )
1149          surf_lsm_h%gtsws = 0.0_wp
1150       ENDIF
1151!--    Horizontal surfaces: urban type         
1152       IF ( urban_surface )  THEN
1153          ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngast ) )
1154          surf_usm_h%gtsws = 0.0_wp
1155       ENDIF
1156!
1157!--    Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
1158!--    westward (l=3) facing
1159       DO  l = 0, 3     
1160          ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngast ) )
1161          surf_def_v(l)%gtsws = 0.0_wp
1162          IF ( land_surface )  THEN
1163             ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngast ) )
1164             surf_lsm_v(l)%gtsws = 0.0_wp
1165          ENDIF
1166          IF ( urban_surface )  THEN
1167             ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngast ) )
1168             surf_usm_v(l)%gtsws = 0.0_wp
1169          ENDIF
1170       ENDDO
1171    ENDIF
1172   
1173 END SUBROUTINE salsa_init_arrays
1174
1175!------------------------------------------------------------------------------!
1176! Description:
1177! ------------
1178!> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA.
1179!> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are
1180!> also merged here.
1181!------------------------------------------------------------------------------!
1182 SUBROUTINE salsa_init
1183
1184    IMPLICIT NONE
1185   
1186    INTEGER(iwp) :: b
1187    INTEGER(iwp) :: c
1188    INTEGER(iwp) :: g
1189    INTEGER(iwp) :: i
1190    INTEGER(iwp) :: j
1191   
1192    CALL location_message( 'initializing SALSA model', .TRUE. )
1193   
1194    bin_low_limits = 0.0_wp
1195    nsect          = 0.0_wp
1196    massacc        = 1.0_wp 
1197   
1198!
1199!-- Indices for chemical components used (-1 = not used)
1200    i = 0
1201    IF ( is_used( prtcl, 'SO4' ) )  THEN
1202       iso4 = get_index( prtcl,'SO4' )
1203       i = i + 1
1204    ENDIF
1205    IF ( is_used( prtcl,'OC' ) )  THEN
1206       ioc = get_index(prtcl, 'OC')
1207       i = i + 1
1208    ENDIF
1209    IF ( is_used( prtcl, 'BC' ) )  THEN
1210       ibc = get_index( prtcl, 'BC' )
1211       i = i + 1
1212    ENDIF
1213    IF ( is_used( prtcl, 'DU' ) )  THEN
1214       idu = get_index( prtcl, 'DU' )
1215       i = i + 1
1216    ENDIF
1217    IF ( is_used( prtcl, 'SS' ) )  THEN
1218       iss = get_index( prtcl, 'SS' )
1219       i = i + 1
1220    ENDIF
1221    IF ( is_used( prtcl, 'NO' ) )  THEN
1222       ino = get_index( prtcl, 'NO' )
1223       i = i + 1
1224    ENDIF
1225    IF ( is_used( prtcl, 'NH' ) )  THEN
1226       inh = get_index( prtcl, 'NH' )
1227       i = i + 1
1228    ENDIF
1229!   
1230!-- All species must be known
1231    IF ( i /= ncc )  THEN
1232       message_string = 'Unknown aerosol species/component(s) given in the' // &
1233                        ' initialization'
1234       CALL message( 'salsa_mod: salsa_init', 'SA0020', 1, 2, 0, 6, 0 )
1235    ENDIF
1236   
1237!
1238!-- Initialise
1239!
1240!-- Aerosol size distribution (TYPE t_section)
1241    aero(:)%dwet     = 1.0E-10_wp
1242    aero(:)%veqh2o   = 1.0E-10_wp
1243    aero(:)%numc     = nclim
1244    aero(:)%core     = 1.0E-10_wp
1245    DO c = 1, maxspec+1    ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
1246       aero(:)%volc(c) = 0.0_wp
1247    ENDDO
1248   
1249    IF ( nldepo )  sedim_vd = 0.0_wp 
1250   
1251    DO  b = 1, nbins
1252       IF ( .NOT. read_restart_data_salsa )  aerosol_number(b)%conc = nclim
1253       aerosol_number(b)%conc_p    = 0.0_wp
1254       aerosol_number(b)%tconc_m   = 0.0_wp
1255       aerosol_number(b)%flux_s    = 0.0_wp
1256       aerosol_number(b)%diss_s    = 0.0_wp
1257       aerosol_number(b)%flux_l    = 0.0_wp
1258       aerosol_number(b)%diss_l    = 0.0_wp
1259       aerosol_number(b)%init      = nclim
1260       aerosol_number(b)%sums_ws_l = 0.0_wp
1261    ENDDO
1262    DO  c = 1, ncc_tot*nbins
1263       IF ( .NOT. read_restart_data_salsa )  aerosol_mass(c)%conc = mclim
1264       aerosol_mass(c)%conc_p    = 0.0_wp
1265       aerosol_mass(c)%tconc_m   = 0.0_wp
1266       aerosol_mass(c)%flux_s    = 0.0_wp
1267       aerosol_mass(c)%diss_s    = 0.0_wp
1268       aerosol_mass(c)%flux_l    = 0.0_wp
1269       aerosol_mass(c)%diss_l    = 0.0_wp
1270       aerosol_mass(c)%init      = mclim
1271       aerosol_mass(c)%sums_ws_l = 0.0_wp
1272    ENDDO
1273   
1274    IF ( .NOT. salsa_gases_from_chem )  THEN
1275       DO  g = 1, ngast
1276          salsa_gas(g)%conc_p    = 0.0_wp
1277          salsa_gas(g)%tconc_m   = 0.0_wp
1278          salsa_gas(g)%flux_s    = 0.0_wp
1279          salsa_gas(g)%diss_s    = 0.0_wp
1280          salsa_gas(g)%flux_l    = 0.0_wp
1281          salsa_gas(g)%diss_l    = 0.0_wp
1282          salsa_gas(g)%sums_ws_l = 0.0_wp
1283       ENDDO
1284       IF ( .NOT. read_restart_data_salsa )  THEN
1285          salsa_gas(1)%conc = H2SO4_init
1286          salsa_gas(2)%conc = HNO3_init
1287          salsa_gas(3)%conc = NH3_init
1288          salsa_gas(4)%conc = OCNV_init
1289          salsa_gas(5)%conc = OCSV_init 
1290       ENDIF
1291!
1292!--    Set initial value for gas compound tracers and initial values
1293       salsa_gas(1)%init = H2SO4_init
1294       salsa_gas(2)%init = HNO3_init
1295       salsa_gas(3)%init = NH3_init
1296       salsa_gas(4)%init = OCNV_init
1297       salsa_gas(5)%init = OCSV_init     
1298    ENDIF
1299!
1300!-- Aerosol radius in each bin: dry and wet (m)
1301    Ra_dry = 1.0E-10_wp
1302!   
1303!-- Initialise aerosol tracers   
1304    aero(:)%vhilim   = 0.0_wp
1305    aero(:)%vlolim   = 0.0_wp
1306    aero(:)%vratiohi = 0.0_wp
1307    aero(:)%vratiolo = 0.0_wp
1308    aero(:)%dmid     = 0.0_wp
1309!
1310!-- Initialise the sectional particle size distribution
1311    CALL set_sizebins()
1312!
1313!-- Initialise location-dependent aerosol size distributions and
1314!-- chemical compositions:
1315    CALL aerosol_init 
1316!
1317!-- Initalisation run of SALSA
1318    DO  i = nxl, nxr
1319       DO  j = nys, nyn
1320          CALL salsa_driver( i, j, 1 )
1321          CALL salsa_diagnostics( i, j )
1322       ENDDO
1323    ENDDO 
1324!
1325!-- Set the aerosol and gas sources
1326    IF ( salsa_source_mode == 'read_from_file' )  THEN
1327       CALL salsa_set_source
1328    ENDIF
1329   
1330    CALL location_message( 'finished', .TRUE. )
1331   
1332 END SUBROUTINE salsa_init
1333
1334!------------------------------------------------------------------------------!
1335! Description:
1336! ------------
1337!> Initializes particle size distribution grid by calculating size bin limits
1338!> and mid-size for *dry* particles in each bin. Called from salsa_initialize
1339!> (only at the beginning of simulation).
1340!> Size distribution described using:
1341!>   1) moving center method (subranges 1 and 2)
1342!>      (Jacobson, Atmos. Env., 31, 131-144, 1997)
1343!>   2) fixed sectional method (subrange 3)
1344!> Size bins in each subrange are spaced logarithmically
1345!> based on given subrange size limits and bin number.
1346!
1347!> Mona changed 06/2017: Use geometric mean diameter to describe the mean
1348!> particle diameter in a size bin, not the arithmeric mean which clearly
1349!> overestimates the total particle volume concentration.
1350!
1351!> Coded by:
1352!> Hannele Korhonen (FMI) 2005
1353!> Harri Kokkola (FMI) 2006
1354!
1355!> Bug fixes for box model + updated for the new aerosol datatype:
1356!> Juha Tonttila (FMI) 2014
1357!------------------------------------------------------------------------------!
1358 SUBROUTINE set_sizebins
1359               
1360    IMPLICIT NONE
1361!   
1362!-- Local variables
1363    INTEGER(iwp) ::  cc
1364    INTEGER(iwp) ::  dd
1365    REAL(wp) ::  ratio_d !< ratio of the upper and lower diameter of subranges
1366!
1367!-- vlolim&vhilim: min & max *dry* volumes [fxm]
1368!-- dmid: bin mid *dry* diameter (m)
1369!-- vratiolo&vratiohi: volume ratio between the center and low/high limit
1370!
1371!-- 1) Size subrange 1:
1372    ratio_d = reglim(2) / reglim(1)   ! section spacing (m)
1373    DO  cc = in1a,fn1a
1374       aero(cc)%vlolim = api6 * ( reglim(1) * ratio_d **                       &
1375                                ( REAL( cc-1 ) / nbin(1) ) ) ** 3.0_wp
1376       aero(cc)%vhilim = api6 * ( reglim(1) * ratio_d **                       &
1377                                ( REAL( cc ) / nbin(1) ) ) ** 3.0_wp
1378       aero(cc)%dmid = SQRT( ( aero(cc)%vhilim / api6 ) ** ( 1.0_wp / 3.0_wp ) &
1379                           * ( aero(cc)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) )
1380       aero(cc)%vratiohi = aero(cc)%vhilim / ( api6 * aero(cc)%dmid ** 3.0_wp )
1381       aero(cc)%vratiolo = aero(cc)%vlolim / ( api6 * aero(cc)%dmid ** 3.0_wp )
1382    ENDDO
1383!
1384!-- 2) Size subrange 2:
1385!-- 2.1) Sub-subrange 2a: high hygroscopicity
1386    ratio_d = reglim(3) / reglim(2)   ! section spacing
1387    DO  dd = in2a, fn2a
1388       cc = dd - in2a
1389       aero(dd)%vlolim = api6 * ( reglim(2) * ratio_d **                       &
1390                                  ( REAL( cc ) / nbin(2) ) ) ** 3.0_wp
1391       aero(dd)%vhilim = api6 * ( reglim(2) * ratio_d **                       &
1392                                  ( REAL( cc+1 ) / nbin(2) ) ) ** 3.0_wp
1393       aero(dd)%dmid = SQRT( ( aero(dd)%vhilim / api6 ) ** ( 1.0_wp / 3.0_wp ) &
1394                           * ( aero(dd)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) )
1395       aero(dd)%vratiohi = aero(dd)%vhilim / ( api6 * aero(dd)%dmid ** 3.0_wp )
1396       aero(dd)%vratiolo = aero(dd)%vlolim / ( api6 * aero(dd)%dmid ** 3.0_wp )
1397    ENDDO
1398!         
1399!-- 2.2) Sub-subrange 2b: low hygroscopicity
1400    IF ( .NOT. no_insoluble )  THEN
1401       aero(in2b:fn2b)%vlolim   = aero(in2a:fn2a)%vlolim
1402       aero(in2b:fn2b)%vhilim   = aero(in2a:fn2a)%vhilim
1403       aero(in2b:fn2b)%dmid     = aero(in2a:fn2a)%dmid
1404       aero(in2b:fn2b)%vratiohi = aero(in2a:fn2a)%vratiohi
1405       aero(in2b:fn2b)%vratiolo = aero(in2a:fn2a)%vratiolo
1406    ENDIF
1407!         
1408!-- Initialize the wet diameter with the bin dry diameter to avoid numerical
1409!-- problems later
1410    aero(:)%dwet = aero(:)%dmid
1411!
1412!-- Save bin limits (lower diameter) to be delivered to the host model if needed
1413    DO cc = 1, nbins
1414       bin_low_limits(cc) = ( aero(cc)%vlolim / api6 )**( 1.0_wp / 3.0_wp )
1415    ENDDO   
1416   
1417 END SUBROUTINE set_sizebins
1418 
1419!------------------------------------------------------------------------------!
1420! Description:
1421! ------------
1422!> Initilize altitude-dependent aerosol size distributions and compositions.
1423!>
1424!> Mona added 06/2017: Correct the number and mass concentrations by normalizing
1425!< by the given total number and mass concentration.
1426!>
1427!> Tomi Raatikainen, FMI, 29.2.2016
1428!------------------------------------------------------------------------------!
1429 SUBROUTINE aerosol_init
1430 
1431    USE arrays_3d,                                                             &
1432        ONLY:  zu
1433 
1434!    USE NETCDF
1435   
1436    USE netcdf_data_input_mod,                                                 &
1437        ONLY:  get_attribute, get_variable,                                    &
1438               netcdf_data_input_get_dimension_length, open_read_file
1439   
1440    IMPLICIT NONE
1441   
1442    INTEGER(iwp) ::  b          !< loop index: size bins
1443    INTEGER(iwp) ::  c          !< loop index: chemical components
1444    INTEGER(iwp) ::  ee         !< index: end
1445    INTEGER(iwp) ::  g          !< loop index: gases
1446    INTEGER(iwp) ::  i          !< loop index: x-direction
1447    INTEGER(iwp) ::  id_faero   !< NetCDF id of PIDS_SALSA
1448    INTEGER(iwp) ::  id_fchem   !< NetCDF id of PIDS_CHEM
1449    INTEGER(iwp) ::  j          !< loop index: y-direction
1450    INTEGER(iwp) ::  k          !< loop index: z-direction
1451    INTEGER(iwp) ::  kk         !< loop index: z-direction
1452    INTEGER(iwp) ::  nz_file    !< Number of grid-points in file (heights)                           
1453    INTEGER(iwp) ::  prunmode
1454    INTEGER(iwp) ::  ss !< index: start
1455    LOGICAL  ::  netcdf_extend = .FALSE. !< Flag indicating wether netcdf
1456                                         !< topography input file or not
1457    REAL(wp), DIMENSION(nbins) ::  core  !< size of the bin mid aerosol particle,
1458    REAL(wp) ::  flag           !< flag to mask topography grid points
1459    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_gas !< gas profiles
1460    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_a !< mass fraction
1461                                                              !< profiles: a
1462    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_b !< and b
1463    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_nsect !< sectional size
1464                                                       !< distribution profile
1465    REAL(wp), DIMENSION(nbins)            ::  nsect  !< size distribution (#/m3)
1466    REAL(wp), DIMENSION(0:nz+1,nbins)     ::  pndist !< size dist as a function
1467                                                     !< of height (#/m3)
1468    REAL(wp), DIMENSION(0:nz+1)           ::  pnf2a  !< number fraction: bins 2a
1469    REAL(wp), DIMENSION(0:nz+1,maxspec)   ::  pvf2a  !< mass distributions of 
1470                                                     !< aerosol species for a 
1471    REAL(wp), DIMENSION(0:nz+1,maxspec)   ::  pvf2b  !< and b-bins     
1472    REAL(wp), DIMENSION(0:nz+1)           ::  pvfOC1a !< mass fraction between
1473                                                     !< SO4 and OC in 1a
1474    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pr_z
1475
1476    prunmode = 1
1477!
1478!-- Bin mean aerosol particle volume (m3)
1479    core(:) = 0.0_wp
1480    core(1:nbins) = api6 * aero(1:nbins)%dmid ** 3.0_wp
1481!   
1482!-- Set concentrations to zero
1483    nsect(:)     = 0.0_wp
1484    pndist(:,:)  = 0.0_wp
1485    pnf2a(:)     = nf2a   
1486    pvf2a(:,:)   = 0.0_wp
1487    pvf2b(:,:)   = 0.0_wp
1488    pvfOC1a(:)   = 0.0_wp
1489
1490    IF ( isdtyp == 1 )  THEN
1491!
1492!--    Read input profiles from PIDS_SALSA   
1493#if defined( __netcdf )
1494!   
1495!--    Location-dependent size distributions and compositions.     
1496       INQUIRE( FILE='PIDS_SALSA'// TRIM( coupling_char ), EXIST=netcdf_extend )
1497       IF ( netcdf_extend )  THEN
1498!
1499!--       Open file in read-only mode 
1500          CALL open_read_file( 'PIDS_SALSA' // TRIM( coupling_char ), id_faero )
1501!
1502!--       Input heights   
1503          CALL netcdf_data_input_get_dimension_length( id_faero, nz_file,      &
1504                                                       "profile_z" ) 
1505         
1506          ALLOCATE( pr_z(nz_file), pr_mass_fracs_a(maxspec,nz_file),           &
1507                    pr_mass_fracs_b(maxspec,nz_file), pr_nsect(nbins,nz_file) ) 
1508          CALL get_variable( id_faero, 'profile_z', pr_z ) 
1509!       
1510!--       Mass fracs profile: 1: H2SO4 (sulphuric acid), 2: OC (organic carbon),
1511!--                           3: BC (black carbon),      4: DU (dust), 
1512!--                           5: SS (sea salt),          6: HNO3 (nitric acid),
1513!--                           7: NH3 (ammonia)         
1514          CALL get_variable( id_faero, "profile_mass_fracs_a", pr_mass_fracs_a,&
1515                             0, nz_file-1, 0, maxspec-1 )
1516          CALL get_variable( id_faero, "profile_mass_fracs_b", pr_mass_fracs_b,&
1517                             0, nz_file-1, 0, maxspec-1 )
1518          CALL get_variable( id_faero, "profile_nsect", pr_nsect, 0, nz_file-1,&
1519                             0, nbins-1 )                   
1520         
1521          kk = 1
1522          DO  k = nzb, nz+1
1523             IF ( kk < nz_file )  THEN
1524                DO  WHILE ( pr_z(kk+1) <= zu(k) )
1525                   kk = kk + 1
1526                   IF ( kk == nz_file )  EXIT
1527                ENDDO
1528             ENDIF
1529             IF ( kk < nz_file )  THEN
1530!             
1531!--             Set initial value for gas compound tracers and initial values
1532                pvf2a(k,:) = pr_mass_fracs_a(:,kk) + ( zu(k) - pr_z(kk) ) / (  &
1533                            pr_z(kk+1) - pr_z(kk) ) * ( pr_mass_fracs_a(:,kk+1)&
1534                            - pr_mass_fracs_a(:,kk) )   
1535                pvf2b(k,:) = pr_mass_fracs_b(:,kk) + ( zu(k) - pr_z(kk) ) / (  &
1536                            pr_z(kk+1) - pr_z(kk) ) * ( pr_mass_fracs_b(:,kk+1)&
1537                            - pr_mass_fracs_b(:,kk) )             
1538                pndist(k,:) = pr_nsect(:,kk) + ( zu(k) - pr_z(kk) ) / (        &
1539                              pr_z(kk+1) - pr_z(kk) ) * ( pr_nsect(:,kk+1) -   &
1540                              pr_nsect(:,kk) )
1541             ELSE
1542                pvf2a(k,:) = pr_mass_fracs_a(:,kk)       
1543                pvf2b(k,:) = pr_mass_fracs_b(:,kk)
1544                pndist(k,:) = pr_nsect(:,kk)
1545             ENDIF
1546             IF ( iso4 < 0 )  THEN
1547                pvf2a(k,1) = 0.0_wp
1548                pvf2b(k,1) = 0.0_wp
1549             ENDIF
1550             IF ( ioc < 0 )  THEN
1551                pvf2a(k,2) = 0.0_wp
1552                pvf2b(k,2) = 0.0_wp
1553             ENDIF
1554             IF ( ibc < 0 )  THEN
1555                pvf2a(k,3) = 0.0_wp
1556                pvf2b(k,3) = 0.0_wp
1557             ENDIF
1558             IF ( idu < 0 )  THEN
1559                pvf2a(k,4) = 0.0_wp
1560                pvf2b(k,4) = 0.0_wp
1561             ENDIF
1562             IF ( iss < 0 )  THEN
1563                pvf2a(k,5) = 0.0_wp
1564                pvf2b(k,5) = 0.0_wp
1565             ENDIF
1566             IF ( ino < 0 )  THEN
1567                pvf2a(k,6) = 0.0_wp
1568                pvf2b(k,6) = 0.0_wp
1569             ENDIF
1570             IF ( inh < 0 )  THEN
1571                pvf2a(k,7) = 0.0_wp
1572                pvf2b(k,7) = 0.0_wp
1573             ENDIF
1574!
1575!--          Then normalise the mass fraction so that SUM = 1
1576             pvf2a(k,:) = pvf2a(k,:) / SUM( pvf2a(k,:) )
1577             IF ( SUM( pvf2b(k,:) ) > 0.0_wp ) pvf2b(k,:) = pvf2b(k,:) /       &
1578                                                            SUM( pvf2b(k,:) )
1579          ENDDO         
1580          DEALLOCATE( pr_z, pr_mass_fracs_a, pr_mass_fracs_b, pr_nsect )
1581       ELSE
1582          message_string = 'Input file '// TRIM( 'PIDS_SALSA' ) //             &
1583                           TRIM( coupling_char ) // ' for SALSA missing!'
1584          CALL message( 'salsa_mod: aerosol_init', 'SA0032', 1, 2, 0, 6, 0 )               
1585       ENDIF   ! netcdf_extend   
1586#endif
1587 
1588    ELSEIF ( isdtyp == 0 )  THEN
1589!
1590!--    Mass fractions for species in a and b-bins
1591       IF ( iso4 > 0 )  THEN
1592          pvf2a(:,1) = mass_fracs_a(iso4) 
1593          pvf2b(:,1) = mass_fracs_b(iso4)
1594       ENDIF
1595       IF ( ioc > 0 )  THEN
1596          pvf2a(:,2) = mass_fracs_a(ioc)
1597          pvf2b(:,2) = mass_fracs_b(ioc) 
1598       ENDIF
1599       IF ( ibc > 0 )  THEN
1600          pvf2a(:,3) = mass_fracs_a(ibc) 
1601          pvf2b(:,3) = mass_fracs_b(ibc)
1602       ENDIF
1603       IF ( idu > 0 )  THEN
1604          pvf2a(:,4) = mass_fracs_a(idu)
1605          pvf2b(:,4) = mass_fracs_b(idu) 
1606       ENDIF
1607       IF ( iss > 0 )  THEN
1608          pvf2a(:,5) = mass_fracs_a(iss)
1609          pvf2b(:,5) = mass_fracs_b(iss) 
1610       ENDIF
1611       IF ( ino > 0 )  THEN
1612          pvf2a(:,6) = mass_fracs_a(ino)
1613          pvf2b(:,6) = mass_fracs_b(ino)
1614       ENDIF
1615       IF ( inh > 0 )  THEN
1616          pvf2a(:,7) = mass_fracs_a(inh)
1617          pvf2b(:,7) = mass_fracs_b(inh)
1618       ENDIF
1619       DO  k = nzb, nz+1
1620          pvf2a(k,:) = pvf2a(k,:) / SUM( pvf2a(k,:) )
1621          IF ( SUM( pvf2b(k,:) ) > 0.0_wp ) pvf2b(k,:) = pvf2b(k,:) /          &
1622                                                         SUM( pvf2b(k,:) )
1623       ENDDO
1624       
1625       CALL size_distribution( n_lognorm, dpg, sigmag, nsect )
1626!
1627!--    Normalize by the given total number concentration
1628       nsect = nsect * SUM( n_lognorm ) * 1.0E+6_wp / SUM( nsect )     
1629       DO  b = in1a, fn2b
1630          pndist(:,b) = nsect(b)
1631       ENDDO
1632    ENDIF
1633   
1634    IF ( igctyp == 1 )  THEN
1635!
1636!--    Read input profiles from PIDS_CHEM   
1637#if defined( __netcdf )
1638!   
1639!--    Location-dependent size distributions and compositions.     
1640       INQUIRE( FILE='PIDS_CHEM' // TRIM( coupling_char ), EXIST=netcdf_extend )
1641       IF ( netcdf_extend  .AND.  .NOT. salsa_gases_from_chem )  THEN
1642!
1643!--       Open file in read-only mode     
1644          CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_fchem )
1645!
1646!--       Input heights   
1647          CALL netcdf_data_input_get_dimension_length( id_fchem, nz_file,      &
1648                                                       "profile_z" ) 
1649          ALLOCATE( pr_z(nz_file), pr_gas(ngast,nz_file) ) 
1650          CALL get_variable( id_fchem, 'profile_z', pr_z ) 
1651!       
1652!--       Gases:
1653          CALL get_variable( id_fchem, "profile_H2SO4", pr_gas(1,:) )
1654          CALL get_variable( id_fchem, "profile_HNO3", pr_gas(2,:) )
1655          CALL get_variable( id_fchem, "profile_NH3", pr_gas(3,:) )
1656          CALL get_variable( id_fchem, "profile_OCNV", pr_gas(4,:) )
1657          CALL get_variable( id_fchem, "profile_OCSV", pr_gas(5,:) )
1658         
1659          kk = 1
1660          DO  k = nzb, nz+1
1661             IF ( kk < nz_file )  THEN
1662                DO  WHILE ( pr_z(kk+1) <= zu(k) )
1663                   kk = kk + 1
1664                   IF ( kk == nz_file )  EXIT
1665                ENDDO
1666             ENDIF
1667             IF ( kk < nz_file )  THEN
1668!             
1669!--             Set initial value for gas compound tracers and initial values
1670                DO  g = 1, ngast
1671                   salsa_gas(g)%init(k) =  pr_gas(g,kk) + ( zu(k) - pr_z(kk) ) &
1672                                           / ( pr_z(kk+1) - pr_z(kk) ) *       &
1673                                           ( pr_gas(g,kk+1) - pr_gas(g,kk) )
1674                   salsa_gas(g)%conc(k,:,:) = salsa_gas(g)%init(k)
1675                ENDDO
1676             ELSE
1677                DO  g = 1, ngast
1678                   salsa_gas(g)%init(k) =  pr_gas(g,kk) 
1679                   salsa_gas(g)%conc(k,:,:) = salsa_gas(g)%init(k)
1680                ENDDO
1681             ENDIF
1682          ENDDO
1683         
1684          DEALLOCATE( pr_z, pr_gas )
1685       ELSEIF ( .NOT. netcdf_extend  .AND.  .NOT.  salsa_gases_from_chem )  THEN
1686          message_string = 'Input file '// TRIM( 'PIDS_CHEM' ) //              &
1687                           TRIM( coupling_char ) // ' for SALSA missing!'
1688          CALL message( 'salsa_mod: aerosol_init', 'SA0033', 1, 2, 0, 6, 0 )               
1689       ENDIF   ! netcdf_extend     
1690#endif
1691
1692    ENDIF
1693
1694    IF ( ioc > 0  .AND.  iso4 > 0 )  THEN     
1695!--    Both are there, so use the given "massDistrA"
1696       pvfOC1a(:) = pvf2a(:,2) / ( pvf2a(:,2) + pvf2a(:,1) )  ! Normalize
1697    ELSEIF ( ioc > 0 )  THEN
1698!--    Pure organic carbon
1699       pvfOC1a(:) = 1.0_wp
1700    ELSEIF ( iso4 > 0 )  THEN
1701!--    Pure SO4
1702       pvfOC1a(:) = 0.0_wp   
1703    ELSE
1704       message_string = 'Either OC or SO4 must be active for aerosol region 1a!'
1705       CALL message( 'salsa_mod: aerosol_init', 'SA0021', 1, 2, 0, 6, 0 )
1706    ENDIF   
1707   
1708!
1709!-- Initialize concentrations
1710    DO  i = nxlg, nxrg 
1711       DO  j = nysg, nyng
1712          DO  k = nzb, nzt+1
1713!
1714!--          Predetermine flag to mask topography         
1715             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1716!         
1717!--          a) Number concentrations
1718!--           Region 1:
1719             DO  b = in1a, fn1a
1720                aerosol_number(b)%conc(k,j,i) = pndist(k,b) * flag
1721                IF ( prunmode == 1 )  THEN
1722                   aerosol_number(b)%init = pndist(:,b)
1723                ENDIF
1724             ENDDO
1725!             
1726!--           Region 2:
1727             IF ( nreg > 1 )  THEN
1728                DO  b = in2a, fn2a
1729                   aerosol_number(b)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) *   &
1730                                                    pndist(k,b) * flag
1731                   IF ( prunmode == 1 )  THEN
1732                      aerosol_number(b)%init = MAX( 0.0_wp, nf2a ) * pndist(:,b)
1733                   ENDIF
1734                ENDDO
1735                IF ( .NOT. no_insoluble )  THEN
1736                   DO  b = in2b, fn2b
1737                      IF ( pnf2a(k) < 1.0_wp )  THEN             
1738                         aerosol_number(b)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp   &
1739                                               - pnf2a(k) ) * pndist(k,b) * flag
1740                         IF ( prunmode == 1 )  THEN
1741                            aerosol_number(b)%init = MAX( 0.0_wp, 1.0_wp -     &
1742                                                          nf2a ) * pndist(:,b)
1743                         ENDIF
1744                      ENDIF
1745                   ENDDO
1746                ENDIF
1747             ENDIF
1748!
1749!--          b) Aerosol mass concentrations
1750!--             bin subrange 1: done here separately due to the SO4/OC convention
1751!--          SO4:
1752             IF ( iso4 > 0 )  THEN
1753                ss = ( iso4 - 1 ) * nbins + in1a !< start
1754                ee = ( iso4 - 1 ) * nbins + fn1a !< end
1755                b = in1a
1756                DO  c = ss, ee
1757                   aerosol_mass(c)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp -         &
1758                                                  pvfOC1a(k) ) * pndist(k,b) * &
1759                                                  core(b) * arhoh2so4 * flag
1760                   IF ( prunmode == 1 )  THEN
1761                      aerosol_mass(c)%init = MAX( 0.0_wp, 1.0_wp - MAXVAL(     &
1762                                             pvfOC1a ) ) * pndist(:,b) *       &
1763                                             core(b) * arhoh2so4
1764                   ENDIF
1765                   b = b+1
1766                ENDDO
1767             ENDIF
1768!--          OC:
1769             IF ( ioc > 0 ) THEN
1770                ss = ( ioc - 1 ) * nbins + in1a !< start
1771                ee = ( ioc - 1 ) * nbins + fn1a !< end
1772                b = in1a
1773                DO  c = ss, ee 
1774                   aerosol_mass(c)%conc(k,j,i) = MAX( 0.0_wp, pvfOC1a(k) ) *   &
1775                                           pndist(k,b) * core(b) * arhooc * flag
1776                   IF ( prunmode == 1 )  THEN
1777                      aerosol_mass(c)%init = MAX( 0.0_wp, MAXVAL( pvfOC1a ) )  &
1778                                             * pndist(:,b) *  core(b) * arhooc
1779                   ENDIF
1780                   b = b+1
1781                ENDDO 
1782             ENDIF
1783             
1784             prunmode = 3  ! Init only once
1785 
1786          ENDDO !< k
1787       ENDDO !< j
1788    ENDDO !< i
1789   
1790!
1791!-- c) Aerosol mass concentrations
1792!--    bin subrange 2:
1793    IF ( nreg > 1 ) THEN
1794   
1795       IF ( iso4 > 0 ) THEN
1796          CALL set_aero_mass( iso4, pvf2a(:,1), pvf2b(:,1), pnf2a, pndist,     &
1797                              core, arhoh2so4 )
1798       ENDIF
1799       IF ( ioc > 0 ) THEN
1800          CALL set_aero_mass( ioc, pvf2a(:,2), pvf2b(:,2), pnf2a, pndist, core,&
1801                              arhooc )
1802       ENDIF
1803       IF ( ibc > 0 ) THEN
1804          CALL set_aero_mass( ibc, pvf2a(:,3), pvf2b(:,3), pnf2a, pndist, core,&
1805                              arhobc )
1806       ENDIF
1807       IF ( idu > 0 ) THEN
1808          CALL set_aero_mass( idu, pvf2a(:,4), pvf2b(:,4), pnf2a, pndist, core,&
1809                              arhodu )
1810       ENDIF
1811       IF ( iss > 0 ) THEN
1812          CALL set_aero_mass( iss, pvf2a(:,5), pvf2b(:,5), pnf2a, pndist, core,&
1813                              arhoss )
1814       ENDIF
1815       IF ( ino > 0 ) THEN
1816          CALL set_aero_mass( ino, pvf2a(:,6), pvf2b(:,6), pnf2a, pndist, core,&
1817                              arhohno3 )
1818       ENDIF
1819       IF ( inh > 0 ) THEN
1820          CALL set_aero_mass( inh, pvf2a(:,7), pvf2b(:,7), pnf2a, pndist, core,&
1821                              arhonh3 )
1822       ENDIF
1823
1824    ENDIF
1825   
1826 END SUBROUTINE aerosol_init
1827 
1828!------------------------------------------------------------------------------!
1829! Description:
1830! ------------
1831!> Create a lognormal size distribution and discretise to a sectional
1832!> representation.
1833!------------------------------------------------------------------------------!
1834 SUBROUTINE size_distribution( in_ntot, in_dpg, in_sigma, psd_sect )
1835   
1836    IMPLICIT NONE
1837   
1838!-- Log-normal size distribution: modes   
1839    REAL(wp), DIMENSION(:), INTENT(in) ::  in_dpg    !< geometric mean diameter
1840                                                     !< (micrometres)
1841    REAL(wp), DIMENSION(:), INTENT(in) ::  in_ntot   !< number conc. (#/cm3)
1842    REAL(wp), DIMENSION(:), INTENT(in) ::  in_sigma  !< standard deviation
1843    REAL(wp), DIMENSION(:), INTENT(inout) ::  psd_sect !< sectional size
1844                                                       !< distribution
1845    INTEGER(iwp) ::  b          !< running index: bin
1846    INTEGER(iwp) ::  ib         !< running index: iteration
1847    REAL(wp) ::  d1             !< particle diameter (m, dummy)
1848    REAL(wp) ::  d2             !< particle diameter (m, dummy)
1849    REAL(wp) ::  delta_d        !< (d2-d1)/10                                                     
1850    REAL(wp) ::  deltadp        !< bin width
1851    REAL(wp) ::  dmidi          !< ( d1 + d2 ) / 2
1852   
1853    DO  b = in1a, fn2b !< aerosol size bins
1854       psd_sect(b) = 0.0_wp
1855!--    Particle diameter at the low limit (largest in the bin) (m)
1856       d1 = ( aero(b)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp )
1857!--    Particle diameter at the high limit (smallest in the bin) (m)
1858       d2 = ( aero(b)%vhilim / api6 ) ** ( 1.0_wp / 3.0_wp )
1859!--    Span of particle diameter in a bin (m)
1860       delta_d = ( d2 - d1 ) / 10.0_wp
1861!--    Iterate:             
1862       DO  ib = 1, 10
1863          d1 = ( aero(b)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) + ( ib - 1)    &
1864               * delta_d
1865          d2 = d1 + delta_d
1866          dmidi = ( d1 + d2 ) / 2.0_wp
1867          deltadp = LOG10( d2 / d1 )
1868         
1869!--       Size distribution
1870!--       in_ntot = total number, total area, or total volume concentration
1871!--       in_dpg = geometric-mean number, area, or volume diameter
1872!--       n(k) = number, area, or volume concentration in a bin
1873!--       n_lognorm and dpg converted to units of #/m3 and m
1874          psd_sect(b) = psd_sect(b) + SUM( in_ntot * 1.0E+6_wp * deltadp /     &
1875                     ( SQRT( 2.0_wp * pi ) * LOG10( in_sigma ) ) *             &
1876                     EXP( -LOG10( dmidi / ( 1.0E-6_wp * in_dpg ) )**2.0_wp /   &
1877                     ( 2.0_wp * LOG10( in_sigma ) ** 2.0_wp ) ) )
1878 
1879       ENDDO
1880    ENDDO
1881   
1882 END SUBROUTINE size_distribution
1883
1884!------------------------------------------------------------------------------!
1885! Description:
1886! ------------
1887!> Sets the mass concentrations to aerosol arrays in 2a and 2b.
1888!>
1889!> Tomi Raatikainen, FMI, 29.2.2016
1890!------------------------------------------------------------------------------!
1891 SUBROUTINE set_aero_mass( ispec, ppvf2a, ppvf2b, ppnf2a, ppndist, pcore, prho )
1892   
1893    IMPLICIT NONE
1894
1895    INTEGER(iwp), INTENT(in) :: ispec  !< Aerosol species index
1896    REAL(wp), INTENT(in) ::  pcore(nbins) !< Aerosol bin mid core volume   
1897    REAL(wp), INTENT(in) ::  ppndist(0:nz+1,nbins) !< Aerosol size distribution
1898    REAL(wp), INTENT(in) ::  ppnf2a(0:nz+1) !< Number fraction for 2a   
1899    REAL(wp), INTENT(in) ::  ppvf2a(0:nz+1) !< Mass distributions for a
1900    REAL(wp), INTENT(in) ::  ppvf2b(0:nz+1) !< and b bins   
1901    REAL(wp), INTENT(in) ::  prho !< Aerosol density
1902    INTEGER(iwp) ::  b  !< loop index
1903    INTEGER(iwp) ::  c  !< loop index       
1904    INTEGER(iwp) ::  ee !< index: end
1905    INTEGER(iwp) ::  i  !< loop index
1906    INTEGER(iwp) ::  j  !< loop index
1907    INTEGER(iwp) ::  k  !< loop index
1908    INTEGER(iwp) ::  prunmode  !< 1 = initialise
1909    INTEGER(iwp) ::  ss !< index: start
1910    REAL(wp) ::  flag   !< flag to mask topography grid points
1911   
1912    prunmode = 1
1913   
1914    DO i = nxlg, nxrg 
1915       DO j = nysg, nyng
1916          DO k = nzb, nzt+1 
1917!
1918!--          Predetermine flag to mask topography
1919             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
1920!             
1921!--          Regime 2a:
1922             ss = ( ispec - 1 ) * nbins + in2a
1923             ee = ( ispec - 1 ) * nbins + fn2a
1924             b = in2a
1925             DO c = ss, ee
1926                aerosol_mass(c)%conc(k,j,i) = MAX( 0.0_wp, ppvf2a(k) ) *       &
1927                               ppnf2a(k) * ppndist(k,b) * pcore(b) * prho * flag
1928                IF ( prunmode == 1 )  THEN
1929                   aerosol_mass(c)%init = MAX( 0.0_wp, MAXVAL( ppvf2a(:) ) ) * &
1930                                          MAXVAL( ppnf2a ) * pcore(b) * prho * &
1931                                          MAXVAL( ppndist(:,b) ) 
1932                ENDIF
1933                b = b+1
1934             ENDDO
1935!--          Regime 2b:
1936             IF ( .NOT. no_insoluble )  THEN
1937                ss = ( ispec - 1 ) * nbins + in2b
1938                ee = ( ispec - 1 ) * nbins + fn2b
1939                b = in2a
1940                DO c = ss, ee
1941                   aerosol_mass(c)%conc(k,j,i) = MAX( 0.0_wp, ppvf2b(k) ) * (  &
1942                                         1.0_wp - ppnf2a(k) ) * ppndist(k,b) * &
1943                                         pcore(b) * prho * flag
1944                   IF ( prunmode == 1 )  THEN
1945                      aerosol_mass(c)%init = MAX( 0.0_wp, MAXVAL( ppvf2b(:) ) )&
1946                                        * ( 1.0_wp - MAXVAL( ppnf2a ) ) *      &
1947                                        MAXVAL( ppndist(:,b) ) * pcore(b) * prho
1948                   ENDIF
1949                   b = b+1
1950                ENDDO
1951             ENDIF
1952             prunmode = 3  ! Init only once
1953          ENDDO
1954       ENDDO
1955    ENDDO
1956 END SUBROUTINE set_aero_mass
1957
1958!------------------------------------------------------------------------------!
1959! Description:
1960! ------------
1961!> Swapping of timelevels
1962!------------------------------------------------------------------------------!
1963 SUBROUTINE salsa_swap_timelevel( mod_count )
1964
1965    IMPLICIT NONE
1966
1967    INTEGER(iwp), INTENT(IN) ::  mod_count  !<
1968    INTEGER(iwp) ::  b  !<   
1969    INTEGER(iwp) ::  c  !<   
1970    INTEGER(iwp) ::  cc !<
1971    INTEGER(iwp) ::  g  !<
1972
1973    IF ( simulated_time >= time_since_reference_point )  THEN
1974
1975    SELECT CASE ( mod_count )
1976
1977       CASE ( 0 )
1978
1979          DO  b = 1, nbins
1980             aerosol_number(b)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   =>        &
1981                nconc_1(:,:,:,b)
1982             aerosol_number(b)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>        &
1983                nconc_2(:,:,:,b)
1984             DO  c = 1, ncc_tot
1985                cc = ( c-1 ) * nbins + b  ! required due to possible Intel18 bug
1986                aerosol_mass(cc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   =>      &
1987                   mconc_1(:,:,:,cc)
1988                aerosol_mass(cc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>      &
1989                   mconc_2(:,:,:,cc)
1990             ENDDO
1991          ENDDO
1992         
1993          IF ( .NOT. salsa_gases_from_chem )  THEN
1994             DO  g = 1, ngast
1995                salsa_gas(g)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   =>          &
1996                   gconc_1(:,:,:,g)
1997                salsa_gas(g)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>          &
1998                   gconc_2(:,:,:,g)
1999             ENDDO
2000          ENDIF
2001
2002       CASE ( 1 )
2003
2004          DO  b = 1, nbins
2005             aerosol_number(b)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   =>        &
2006                nconc_2(:,:,:,b)
2007             aerosol_number(b)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>        &
2008                nconc_1(:,:,:,b)
2009             DO  c = 1, ncc_tot
2010                cc = ( c-1 ) * nbins + b  ! required due to possible Intel18 bug
2011                aerosol_mass(cc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   =>      &
2012                   mconc_2(:,:,:,cc)
2013                aerosol_mass(cc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>      &
2014                   mconc_1(:,:,:,cc)
2015             ENDDO
2016          ENDDO
2017         
2018          IF ( .NOT. salsa_gases_from_chem )  THEN
2019             DO  g = 1, ngast
2020                salsa_gas(g)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   =>          &
2021                   gconc_2(:,:,:,g)
2022                salsa_gas(g)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>          &
2023                   gconc_1(:,:,:,g)
2024             ENDDO
2025          ENDIF
2026
2027    END SELECT
2028
2029    ENDIF
2030
2031 END SUBROUTINE salsa_swap_timelevel
2032
2033
2034!------------------------------------------------------------------------------!
2035! Description:
2036! ------------
2037!> This routine reads the respective restart data.
2038!------------------------------------------------------------------------------!
2039 SUBROUTINE salsa_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,           &
2040                             nxr_on_file, nynf, nync, nyn_on_file, nysf,       &
2041                             nysc, nys_on_file, tmp_3d, found )
2042
2043   
2044    IMPLICIT NONE
2045   
2046    INTEGER(iwp) ::  b  !<   
2047    INTEGER(iwp) ::  c  !<
2048    INTEGER(iwp) ::  g  !<
2049    INTEGER(iwp) ::  k  !<
2050    INTEGER(iwp) ::  nxlc            !<
2051    INTEGER(iwp) ::  nxlf            !<
2052    INTEGER(iwp) ::  nxl_on_file     !<
2053    INTEGER(iwp) ::  nxrc            !<
2054    INTEGER(iwp) ::  nxrf            !<
2055    INTEGER(iwp) ::  nxr_on_file     !<
2056    INTEGER(iwp) ::  nync            !<
2057    INTEGER(iwp) ::  nynf            !<
2058    INTEGER(iwp) ::  nyn_on_file     !<
2059    INTEGER(iwp) ::  nysc            !<
2060    INTEGER(iwp) ::  nysf            !<
2061    INTEGER(iwp) ::  nys_on_file     !<
2062
2063    LOGICAL, INTENT(OUT)  ::  found
2064
2065    REAL(wp), &
2066       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
2067       
2068    found = .FALSE.
2069   
2070    IF ( read_restart_data_salsa )  THEN
2071   
2072       SELECT CASE ( restart_string(1:length) )
2073       
2074          CASE ( 'aerosol_number' )
2075             DO  b = 1, nbins
2076                IF ( k == 1 )  READ ( 13 ) tmp_3d
2077                aerosol_number(b)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 
2078                               tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2079                found = .TRUE.
2080             ENDDO
2081       
2082          CASE ( 'aerosol_mass' )
2083             DO  c = 1, ncc_tot * nbins
2084                IF ( k == 1 )  READ ( 13 ) tmp_3d
2085                aerosol_mass(c)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 
2086                               tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2087                found = .TRUE.
2088             ENDDO
2089         
2090          CASE ( 'salsa_gas' )
2091             DO  g = 1, ngast
2092                IF ( k == 1 )  READ ( 13 ) tmp_3d
2093                salsa_gas(g)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  & 
2094                               tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2095                found = .TRUE.
2096             ENDDO
2097             
2098          CASE DEFAULT
2099             found = .FALSE.
2100             
2101       END SELECT
2102    ENDIF
2103
2104 END SUBROUTINE salsa_rrd_local
2105   
2106
2107!------------------------------------------------------------------------------!
2108! Description:
2109! ------------
2110!> This routine writes the respective restart data.
2111!> Note that the following input variables in PARIN have to be equal between
2112!> restart runs:
2113!>    listspec, nbin, nbin2, nf2a, ncc, mass_fracs_a, mass_fracs_b
2114!------------------------------------------------------------------------------!
2115 SUBROUTINE salsa_wrd_local
2116
2117    IMPLICIT NONE
2118   
2119    INTEGER(iwp) ::  b  !<   
2120    INTEGER(iwp) ::  c  !<
2121    INTEGER(iwp) ::  g  !<
2122   
2123    IF ( write_binary  .AND.  write_binary_salsa )  THEN
2124   
2125       CALL wrd_write_string( 'aerosol_number' )
2126       DO  b = 1, nbins
2127          WRITE ( 14 )  aerosol_number(b)%conc
2128       ENDDO
2129       
2130       CALL wrd_write_string( 'aerosol_mass' )
2131       DO  c = 1, nbins*ncc_tot
2132          WRITE ( 14 )  aerosol_mass(c)%conc
2133       ENDDO
2134       
2135       CALL wrd_write_string( 'salsa_gas' )
2136       DO  g = 1, ngast
2137          WRITE ( 14 )  salsa_gas(g)%conc
2138       ENDDO
2139         
2140    ENDIF
2141       
2142 END SUBROUTINE salsa_wrd_local   
2143
2144
2145!------------------------------------------------------------------------------!
2146! Description:
2147! ------------
2148!> Performs necessary unit and dimension conversion between the host model and
2149!> SALSA module, and calls the main SALSA routine.
2150!> Partially adobted form the original SALSA boxmodel version.
2151!> Now takes masses in as kg/kg from LES!! Converted to m3/m3 for SALSA
2152!> 05/2016 Juha: This routine is still pretty much in its original shape.
2153!>               It's dumb as a mule and twice as ugly, so implementation of
2154!>               an improved solution is necessary sooner or later.
2155!> Juha Tonttila, FMI, 2014
2156!> Jaakko Ahola, FMI, 2016
2157!> Only aerosol processes included, Mona Kurppa, UHel, 2017
2158!------------------------------------------------------------------------------!
2159 SUBROUTINE salsa_driver( i, j, prunmode )
2160
2161    USE arrays_3d,                                                             &
2162        ONLY: pt_p, q_p, rho_air_zw, u, v, w
2163       
2164    USE plant_canopy_model_mod,                                                &
2165        ONLY: lad_s
2166       
2167    USE surface_mod,                                                           &
2168        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
2169               surf_usm_v
2170 
2171    IMPLICIT NONE
2172   
2173    INTEGER(iwp), INTENT(in) ::  i   !< loop index
2174    INTEGER(iwp), INTENT(in) ::  j   !< loop index
2175    INTEGER(iwp), INTENT(in) ::  prunmode !< 1: Initialization call
2176                                          !< 2: Spinup period call
2177                                          !< 3: Regular runtime call
2178!-- Local variables
2179    TYPE(t_section), DIMENSION(fn2b) ::  aero_old !< helper array
2180    INTEGER(iwp) ::  bb     !< loop index
2181    INTEGER(iwp) ::  cc     !< loop index
2182    INTEGER(iwp) ::  endi   !< end index
2183    INTEGER(iwp) ::  k_wall !< vertical index of topography top
2184    INTEGER(iwp) ::  k      !< loop index
2185    INTEGER(iwp) ::  l      !< loop index
2186    INTEGER(iwp) ::  nc_h2o !< index of H2O in the prtcl index table
2187    INTEGER(iwp) ::  ss     !< loop index
2188    INTEGER(iwp) ::  str    !< start index
2189    INTEGER(iwp) ::  vc     !< default index in prtcl
2190    REAL(wp) ::  cw_old     !< previous H2O mixing ratio
2191    REAL(wp) ::  flag       !< flag to mask topography grid points
2192    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_adn !< air density (kg/m3)   
2193    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cs  !< H2O sat. vapour conc.
2194    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_cw  !< H2O vapour concentration
2195    REAL(wp) ::  in_lad                       !< leaf area density (m2/m3)
2196    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_p   !< pressure (Pa)     
2197    REAL(wp) ::  in_rh                        !< relative humidity                     
2198    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_t   !< temperature (K)
2199    REAL(wp), DIMENSION(nzb:nzt+1) ::  in_u   !< wind magnitude (m/s)
2200    REAL(wp), DIMENSION(nzb:nzt+1) ::  kvis   !< kinematic viscosity of air(m2/s)                                           
2201    REAL(wp), DIMENSION(nzb:nzt+1,fn2b) ::  Sc      !< particle Schmidt number   
2202    REAL(wp), DIMENSION(nzb:nzt+1,fn2b) ::  vd      !< particle fall seed (m/s,
2203                                                    !< sedimentation velocity)
2204    REAL(wp), DIMENSION(nzb:nzt+1) ::  ppm_to_nconc !< Conversion factor
2205                                                    !< from ppm to #/m3                                                     
2206    REAL(wp) ::  zgso4  !< SO4
2207    REAL(wp) ::  zghno3 !< HNO3
2208    REAL(wp) ::  zgnh3  !< NH3
2209    REAL(wp) ::  zgocnv !< non-volatile OC
2210    REAL(wp) ::  zgocsv !< semi-volatile OC
2211   
2212    aero_old(:)%numc = 0.0_wp
2213    in_adn           = 0.0_wp   
2214    in_cs            = 0.0_wp
2215    in_cw            = 0.0_wp 
2216    in_lad           = 0.0_wp
2217    in_rh            = 0.0_wp
2218    in_p             = 0.0_wp 
2219    in_t             = 0.0_wp 
2220    in_u             = 0.0_wp
2221    kvis             = 0.0_wp
2222    Sc               = 0.0_wp
2223    vd               = 0.0_wp
2224    ppm_to_nconc     = 1.0_wp
2225    zgso4            = nclim
2226    zghno3           = nclim
2227    zgnh3            = nclim
2228    zgocnv           = nclim
2229    zgocsv           = nclim
2230   
2231!       
2232!-- Aerosol number is always set, but mass can be uninitialized
2233    DO cc = 1, nbins
2234       aero(cc)%volc     = 0.0_wp
2235       aero_old(cc)%volc = 0.0_wp
2236    ENDDO 
2237!   
2238!-- Set the salsa runtime config (How to make this more efficient?)
2239    CALL set_salsa_runtime( prunmode )
2240!             
2241!-- Calculate thermodynamic quantities needed in SALSA
2242    CALL salsa_thrm_ij( i, j, p_ij=in_p, temp_ij=in_t, cw_ij=in_cw,            &
2243                        cs_ij=in_cs, adn_ij=in_adn )
2244!
2245!-- Magnitude of wind: needed for deposition
2246    IF ( lsdepo )  THEN
2247       in_u(nzb+1:nzt) = SQRT(                                                 &
2248                   ( 0.5_wp * ( u(nzb+1:nzt,j,i) + u(nzb+1:nzt,j,i+1) ) )**2 + & 
2249                   ( 0.5_wp * ( v(nzb+1:nzt,j,i) + v(nzb+1:nzt,j+1,i) ) )**2 + &
2250                   ( 0.5_wp * ( w(nzb:nzt-1,j,i) + w(nzb+1:nzt,j,  i) ) )**2 )
2251    ENDIF
2252!
2253!-- Calculate conversion factors for gas concentrations
2254    ppm_to_nconc = for_ppm_to_nconc * in_p / in_t
2255!
2256!-- Determine topography-top index on scalar grid
2257    k_wall = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,j,i), 12 ) ),          &
2258                     DIM = 1 ) - 1     
2259               
2260    DO k = nzb+1, nzt
2261!
2262!--    Predetermine flag to mask topography
2263       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2264!       
2265!--    Do not run inside buildings       
2266       IF ( flag == 0.0_wp )  CYCLE   
2267!
2268!--    Wind velocity for dry depositon on vegetation   
2269       IF ( lsdepo_vege  .AND.  plant_canopy  )  THEN
2270          in_lad = lad_s(k-k_wall,j,i)
2271       ENDIF       
2272!
2273!--    For initialization and spinup, limit the RH with the parameter rhlim
2274       IF ( prunmode < 3 ) THEN
2275          in_cw(k) = MIN( in_cw(k), in_cs(k) * rhlim )
2276       ELSE
2277          in_cw(k) = in_cw(k)
2278       ENDIF
2279       cw_old = in_cw(k) !* in_adn(k)
2280!               
2281!--    Set volume concentrations:
2282!--    Sulphate (SO4) or sulphuric acid H2SO4
2283       IF ( iso4 > 0 )  THEN
2284          vc = 1
2285          str = ( iso4-1 ) * nbins + 1    ! start index
2286          endi = iso4 * nbins             ! end index
2287          cc = 1
2288          DO ss = str, endi
2289             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2so4
2290             cc = cc+1
2291          ENDDO
2292          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2293       ENDIF
2294       
2295!--    Organic carbon (OC) compounds
2296       IF ( ioc > 0 )  THEN
2297          vc = 2
2298          str = ( ioc-1 ) * nbins + 1
2299          endi = ioc * nbins
2300          cc = 1
2301          DO ss = str, endi
2302             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhooc 
2303             cc = cc+1
2304          ENDDO
2305          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2306       ENDIF
2307       
2308!--    Black carbon (BC)
2309       IF ( ibc > 0 )  THEN
2310          vc = 3
2311          str = ( ibc-1 ) * nbins + 1 + fn1a
2312          endi = ibc * nbins
2313          cc = 1 + fn1a
2314          DO ss = str, endi
2315             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhobc 
2316             cc = cc+1
2317          ENDDO                   
2318          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2319       ENDIF
2320
2321!--    Dust (DU)
2322       IF ( idu > 0 )  THEN
2323          vc = 4
2324          str = ( idu-1 ) * nbins + 1 + fn1a
2325          endi = idu * nbins
2326          cc = 1 + fn1a
2327          DO ss = str, endi
2328             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhodu 
2329             cc = cc+1
2330          ENDDO
2331          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2332       ENDIF
2333
2334!--    Sea salt (SS)
2335       IF ( iss > 0 )  THEN
2336          vc = 5
2337          str = ( iss-1 ) * nbins + 1 + fn1a
2338          endi = iss * nbins
2339          cc = 1 + fn1a
2340          DO ss = str, endi
2341             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoss 
2342             cc = cc+1
2343          ENDDO
2344          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2345       ENDIF
2346
2347!--    Nitrate (NO(3-)) or nitric acid HNO3
2348       IF ( ino > 0 )  THEN
2349          vc = 6
2350          str = ( ino-1 ) * nbins + 1 
2351          endi = ino * nbins
2352          cc = 1
2353          DO ss = str, endi
2354             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhohno3 
2355             cc = cc+1
2356          ENDDO
2357          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2358       ENDIF
2359
2360!--    Ammonium (NH(4+)) or ammonia NH3
2361       IF ( inh > 0 )  THEN
2362          vc = 7
2363          str = ( inh-1 ) * nbins + 1
2364          endi = inh * nbins
2365          cc = 1
2366          DO ss = str, endi
2367             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhonh3 
2368             cc = cc+1
2369          ENDDO
2370          aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2371       ENDIF
2372
2373!--    Water (always used)
2374       nc_h2o = get_index( prtcl,'H2O' )
2375       vc = 8
2376       str = ( nc_h2o-1 ) * nbins + 1
2377       endi = nc_h2o * nbins
2378       cc = 1
2379       IF ( advect_particle_water )  THEN
2380          DO ss = str, endi
2381             aero(cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2o 
2382             cc = cc+1
2383          ENDDO
2384       ELSE
2385         aero(1:nbins)%volc(vc) = mclim 
2386       ENDIF
2387       aero_old(1:nbins)%volc(vc) = aero(1:nbins)%volc(vc)
2388!
2389!--    Number concentrations (numc) and particle sizes
2390!--    (dwet = wet diameter, core = dry volume)
2391       DO  bb = 1, nbins
2392          aero(bb)%numc = aerosol_number(bb)%conc(k,j,i) 
2393          aero_old(bb)%numc = aero(bb)%numc
2394          IF ( aero(bb)%numc > nclim )  THEN
2395             aero(bb)%dwet = ( SUM( aero(bb)%volc(:) ) / aero(bb)%numc / api6 )&
2396                                **( 1.0_wp / 3.0_wp )
2397             aero(bb)%core = SUM( aero(bb)%volc(1:7) ) / aero(bb)%numc 
2398          ELSE
2399             aero(bb)%dwet = aero(bb)%dmid
2400             aero(bb)%core = api6 * ( aero(bb)%dwet ) ** 3.0_wp
2401          ENDIF
2402       ENDDO
2403!       
2404!--    On EACH call of salsa_driver, calculate the ambient sizes of
2405!--    particles by equilibrating soluble fraction of particles with water
2406!--    using the ZSR method.
2407       in_rh = in_cw(k) / in_cs(k)
2408       IF ( prunmode==1  .OR.  .NOT. advect_particle_water )  THEN
2409          CALL equilibration( in_rh, in_t(k), aero, .TRUE. )
2410       ENDIF
2411!
2412!--    Gaseous tracer concentrations in #/m3
2413       IF ( salsa_gases_from_chem )  THEN       
2414!       
2415!--       Convert concentrations in ppm to #/m3
2416          zgso4  = chem_species(gas_index_chem(1))%conc(k,j,i) * ppm_to_nconc(k)
2417          zghno3 = chem_species(gas_index_chem(2))%conc(k,j,i) * ppm_to_nconc(k)
2418          zgnh3  = chem_species(gas_index_chem(3))%conc(k,j,i) * ppm_to_nconc(k)
2419          zgocnv = chem_species(gas_index_chem(4))%conc(k,j,i) * ppm_to_nconc(k)     
2420          zgocsv = chem_species(gas_index_chem(5))%conc(k,j,i) * ppm_to_nconc(k)                 
2421       ELSE
2422          zgso4  = salsa_gas(1)%conc(k,j,i) 
2423          zghno3 = salsa_gas(2)%conc(k,j,i) 
2424          zgnh3  = salsa_gas(3)%conc(k,j,i) 
2425          zgocnv = salsa_gas(4)%conc(k,j,i) 
2426          zgocsv = salsa_gas(5)%conc(k,j,i)
2427       ENDIF   
2428!
2429!--    ***************************************!
2430!--                   Run SALSA               !
2431!--    ***************************************!
2432       CALL run_salsa( in_p(k), in_cw(k), in_cs(k), in_t(k), in_u(k),          &
2433                       in_adn(k), in_lad, zgso4, zgocnv, zgocsv, zghno3, zgnh3,&
2434                       aero, prtcl, kvis(k), Sc(k,:), vd(k,:), dt_salsa )
2435!--    ***************************************!
2436       IF ( lsdepo ) sedim_vd(k,j,i,:) = vd(k,:)
2437!                           
2438!--    Calculate changes in concentrations
2439       DO bb = 1, nbins
2440          aerosol_number(bb)%conc(k,j,i) = aerosol_number(bb)%conc(k,j,i)      &
2441                                 +  ( aero(bb)%numc - aero_old(bb)%numc ) * flag
2442       ENDDO
2443       
2444       IF ( iso4 > 0 )  THEN
2445          vc = 1
2446          str = ( iso4-1 ) * nbins + 1
2447          endi = iso4 * nbins
2448          cc = 1
2449          DO ss = str, endi
2450             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2451                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2452                               * arhoh2so4 * flag
2453             cc = cc+1
2454          ENDDO
2455       ENDIF
2456       
2457       IF ( ioc > 0 )  THEN
2458          vc = 2
2459          str = ( ioc-1 ) * nbins + 1
2460          endi = ioc * nbins
2461          cc = 1
2462          DO ss = str, endi
2463             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2464                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2465                               * arhooc * flag
2466             cc = cc+1
2467          ENDDO
2468       ENDIF
2469       
2470       IF ( ibc > 0 )  THEN
2471          vc = 3
2472          str = ( ibc-1 ) * nbins + 1 + fn1a
2473          endi = ibc * nbins
2474          cc = 1 + fn1a
2475          DO ss = str, endi
2476             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2477                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2478                               * arhobc * flag 
2479             cc = cc+1
2480          ENDDO
2481       ENDIF
2482       
2483       IF ( idu > 0 )  THEN
2484          vc = 4
2485          str = ( idu-1 ) * nbins + 1 + fn1a
2486          endi = idu * nbins
2487          cc = 1 + fn1a
2488          DO ss = str, endi
2489             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2490                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2491                               * arhodu * flag
2492             cc = cc+1
2493          ENDDO
2494       ENDIF
2495       
2496       IF ( iss > 0 )  THEN
2497          vc = 5
2498          str = ( iss-1 ) * nbins + 1 + fn1a
2499          endi = iss * nbins
2500          cc = 1 + fn1a
2501          DO ss = str, endi
2502             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2503                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2504                               * arhoss * flag
2505             cc = cc+1
2506          ENDDO
2507       ENDIF
2508       
2509       IF ( ino > 0 )  THEN
2510          vc = 6
2511          str = ( ino-1 ) * nbins + 1
2512          endi = ino * nbins
2513          cc = 1
2514          DO ss = str, endi
2515             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2516                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2517                               * arhohno3 * flag
2518             cc = cc+1
2519          ENDDO
2520       ENDIF
2521       
2522       IF ( inh > 0 )  THEN
2523          vc = 7
2524          str = ( ino-1 ) * nbins + 1
2525          endi = ino * nbins
2526          cc = 1
2527          DO ss = str, endi
2528             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2529                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2530                               * arhonh3 * flag
2531             cc = cc+1
2532          ENDDO
2533       ENDIF
2534       
2535       IF ( advect_particle_water )  THEN
2536          nc_h2o = get_index( prtcl,'H2O' )
2537          vc = 8
2538          str = ( nc_h2o-1 ) * nbins + 1
2539          endi = nc_h2o * nbins
2540          cc = 1
2541          DO ss = str, endi
2542             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i)       &
2543                               + ( aero(cc)%volc(vc) - aero_old(cc)%volc(vc) ) &
2544                               * arhoh2o * flag
2545             IF ( prunmode == 1 )  THEN
2546                aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k),      &
2547                                               aerosol_mass(ss)%conc(k,j,i) )
2548             ENDIF
2549             cc = cc+1                             
2550          ENDDO
2551       ENDIF
2552
2553!--    Condensation of precursor gases
2554       IF ( lscndgas )  THEN
2555          IF ( salsa_gases_from_chem )  THEN         
2556!         
2557!--          SO4 (or H2SO4)
2558             chem_species( gas_index_chem(1) )%conc(k,j,i) =                &
2559                            chem_species( gas_index_chem(1) )%conc(k,j,i) + &
2560                                                  ( zgso4 / ppm_to_nconc(k) - &
2561                       chem_species( gas_index_chem(1) )%conc(k,j,i) ) * flag
2562!                           
2563!--          HNO3
2564             chem_species( gas_index_chem(2) )%conc(k,j,i) =                &
2565                            chem_species( gas_index_chem(2) )%conc(k,j,i) + &
2566                                                 ( zghno3 / ppm_to_nconc(k) - &
2567                       chem_species( gas_index_chem(2) )%conc(k,j,i) ) * flag
2568!                           
2569!--          NH3
2570             chem_species( gas_index_chem(3) )%conc(k,j,i) =                &
2571                            chem_species( gas_index_chem(3) )%conc(k,j,i) + &
2572                                                  ( zgnh3 / ppm_to_nconc(k) - &
2573                       chem_species( gas_index_chem(3) )%conc(k,j,i) ) * flag
2574!                           
2575!--          non-volatile OC
2576             chem_species( gas_index_chem(4) )%conc(k,j,i) =                &
2577                            chem_species( gas_index_chem(4) )%conc(k,j,i) + &
2578                                                 ( zgocnv / ppm_to_nconc(k) - &
2579                       chem_species( gas_index_chem(4) )%conc(k,j,i) ) * flag
2580!                           
2581!--          semi-volatile OC
2582             chem_species( gas_index_chem(5) )%conc(k,j,i) =                &
2583                            chem_species( gas_index_chem(5) )%conc(k,j,i) + &
2584                                                 ( zgocsv / ppm_to_nconc(k) - &
2585                       chem_species( gas_index_chem(5) )%conc(k,j,i) ) * flag                 
2586         
2587          ELSE
2588!         
2589!--          SO4 (or H2SO4)
2590             salsa_gas(1)%conc(k,j,i) = salsa_gas(1)%conc(k,j,i) + ( zgso4 -   &
2591                                          salsa_gas(1)%conc(k,j,i) ) * flag
2592!                           
2593!--          HNO3
2594             salsa_gas(2)%conc(k,j,i) = salsa_gas(2)%conc(k,j,i) + ( zghno3 -  &
2595                                          salsa_gas(2)%conc(k,j,i) ) * flag
2596!                           
2597!--          NH3
2598             salsa_gas(3)%conc(k,j,i) = salsa_gas(3)%conc(k,j,i) + ( zgnh3 -   &
2599                                          salsa_gas(3)%conc(k,j,i) ) * flag
2600!                           
2601!--          non-volatile OC
2602             salsa_gas(4)%conc(k,j,i) = salsa_gas(4)%conc(k,j,i) + ( zgocnv -  &
2603                                          salsa_gas(4)%conc(k,j,i) ) * flag
2604!                           
2605!--          semi-volatile OC
2606             salsa_gas(5)%conc(k,j,i) = salsa_gas(5)%conc(k,j,i) + ( zgocsv -  &
2607                                          salsa_gas(5)%conc(k,j,i) ) * flag
2608          ENDIF
2609       ENDIF
2610!               
2611!--    Tendency of water vapour mixing ratio is obtained from the
2612!--    change in RH during SALSA run. This releases heat and changes pt.
2613!--    Assumes no temperature change during SALSA run.
2614!--    q = r / (1+r), Euler method for integration
2615!
2616       IF ( feedback_to_palm )  THEN
2617          q_p(k,j,i) = q_p(k,j,i) + 1.0_wp / ( in_cw(k) * in_adn(k) + 1.0_wp ) &
2618                       ** 2.0_wp * ( in_cw(k) - cw_old ) * in_adn(k) 
2619          pt_p(k,j,i) = pt_p(k,j,i) + alv / c_p * ( in_cw(k) - cw_old ) *      &
2620                        in_adn(k) / ( in_cw(k) / in_adn(k) + 1.0_wp ) ** 2.0_wp&
2621                        * pt_p(k,j,i) / in_t(k)
2622       ENDIF
2623                         
2624    ENDDO   ! k
2625!   
2626!-- Set surfaces and wall fluxes due to deposition 
2627    IF ( lsdepo_topo  .AND.  prunmode == 3 )  THEN
2628       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
2629          CALL depo_topo( i, j, surf_def_h(0), vd, Sc, kvis, in_u, rho_air_zw )
2630          DO  l = 0, 3
2631             CALL depo_topo( i, j, surf_def_v(l), vd, Sc, kvis, in_u,          &
2632                             rho_air_zw**0.0_wp )
2633          ENDDO
2634       ELSE
2635          CALL depo_topo( i, j, surf_usm_h, vd, Sc, kvis, in_u, rho_air_zw )
2636          DO  l = 0, 3
2637             CALL depo_topo( i, j, surf_usm_v(l), vd, Sc, kvis, in_u,          &
2638                             rho_air_zw**0.0_wp )
2639          ENDDO
2640          CALL depo_topo( i, j, surf_lsm_h, vd, Sc, kvis, in_u, rho_air_zw )
2641          DO  l = 0, 3
2642             CALL depo_topo( i, j, surf_lsm_v(l), vd, Sc, kvis, in_u,          &
2643                             rho_air_zw**0.0_wp )
2644          ENDDO
2645       ENDIF
2646    ENDIF
2647   
2648 END SUBROUTINE salsa_driver
2649
2650!------------------------------------------------------------------------------!
2651! Description:
2652! ------------
2653!> The SALSA subroutine
2654!> Modified for the new aerosol datatype,
2655!> Juha Tonttila, FMI, 2014.
2656!> Only aerosol processes included, Mona Kurppa, UHel, 2017
2657!------------------------------------------------------------------------------!   
2658 SUBROUTINE run_salsa( ppres, pcw, pcs, ptemp, mag_u, adn, lad, pc_h2so4,      &
2659                       pc_ocnv, pc_ocsv, pc_hno3, pc_nh3, paero, prtcl, kvis,  &
2660                       Sc, vc, ptstep )
2661
2662    IMPLICIT NONE
2663!
2664!-- Input parameters and variables
2665    REAL(wp), INTENT(in) ::  adn    !< air density (kg/m3)
2666    REAL(wp), INTENT(in) ::  lad    !< leaf area density (m2/m3)
2667    REAL(wp), INTENT(in) ::  mag_u  !< magnitude of wind (m/s)
2668    REAL(wp), INTENT(in) ::  ppres  !< atmospheric pressure at each grid
2669                                    !< point (Pa)
2670    REAL(wp), INTENT(in) ::  ptemp  !< temperature at each grid point (K)
2671    REAL(wp), INTENT(in) ::  ptstep !< time step of salsa processes (s)
2672    TYPE(component_index), INTENT(in) :: prtcl  !< part. component index table
2673!       
2674!-- Input variables that are changed within:
2675    REAL(wp), INTENT(inout) ::  kvis     !< kinematic viscosity of air (m2/s)
2676    REAL(wp), INTENT(inout) ::  Sc(:)    !< particle Schmidt number
2677    REAL(wp), INTENT(inout) ::  vc(:)    !< particle fall speed (m/s,
2678                                         !< sedimentation velocity)
2679!-- Gas phase concentrations at each grid point (#/m3)
2680    REAL(wp), INTENT(inout) ::  pc_h2so4 !< sulphuric acid
2681    REAL(wp), INTENT(inout) ::  pc_hno3  !< nitric acid
2682    REAL(wp), INTENT(inout) ::  pc_nh3   !< ammonia
2683    REAL(wp), INTENT(inout) ::  pc_ocnv  !< nonvolatile OC
2684    REAL(wp), INTENT(inout) ::  pc_ocsv  !< semivolatile OC
2685    REAL(wp), INTENT(inout) ::  pcs      !< Saturation concentration of water
2686                                         !< vapour (kg/m3)
2687    REAL(wp), INTENT(inout) ::  pcw      !< Water vapour concentration (kg/m3)                                                   
2688    TYPE(t_section), INTENT(inout) ::  paero(fn2b) 
2689!
2690!-- Coagulation
2691    IF ( lscoag )   THEN
2692       CALL coagulation( paero, ptstep, ptemp, ppres )
2693    ENDIF
2694!
2695!-- Condensation
2696    IF ( lscnd )   THEN
2697       CALL condensation( paero, pc_h2so4, pc_ocnv, pc_ocsv,  pc_hno3, pc_nh3, &
2698                          pcw, pcs, ptemp, ppres, ptstep, prtcl )
2699    ENDIF   
2700!   
2701!-- Deposition
2702    IF ( lsdepo )  THEN
2703       CALL deposition( paero, ptemp, adn, mag_u, lad, kvis, Sc, vc ) 
2704    ENDIF       
2705!
2706!-- Size distribution bin update
2707!-- Mona: why done 3 times in SALSA-standalone?
2708    IF ( lsdistupdate )   THEN
2709       CALL distr_update( paero )
2710    ENDIF
2711   
2712  END SUBROUTINE run_salsa 
2713 
2714!------------------------------------------------------------------------------!
2715! Description:
2716! ------------
2717!> Set logical switches according to the host model state and user-specified
2718!> NAMELIST options.
2719!> Juha Tonttila, FMI, 2014
2720!> Only aerosol processes included, Mona Kurppa, UHel, 2017
2721!------------------------------------------------------------------------------!
2722 SUBROUTINE set_salsa_runtime( prunmode )
2723 
2724    IMPLICIT NONE
2725   
2726    INTEGER(iwp), INTENT(in) ::  prunmode
2727   
2728    SELECT CASE(prunmode)
2729
2730       CASE(1) !< Initialization
2731          lscoag       = .FALSE.
2732          lscnd        = .FALSE.
2733          lscndgas     = .FALSE.
2734          lscndh2oae   = .FALSE.
2735          lsdepo       = .FALSE.
2736          lsdepo_vege  = .FALSE.
2737          lsdepo_topo  = .FALSE.
2738          lsdistupdate = .TRUE.
2739
2740       CASE(2)  !< Spinup period
2741          lscoag      = ( .FALSE. .AND. nlcoag   )
2742          lscnd       = ( .TRUE.  .AND. nlcnd    )
2743          lscndgas    = ( .TRUE.  .AND. nlcndgas )
2744          lscndh2oae  = ( .TRUE.  .AND. nlcndh2oae )
2745
2746       CASE(3)  !< Run
2747          lscoag       = nlcoag
2748          lscnd        = nlcnd
2749          lscndgas     = nlcndgas
2750          lscndh2oae   = nlcndh2oae
2751          lsdepo       = nldepo
2752          lsdepo_vege  = nldepo_vege
2753          lsdepo_topo  = nldepo_topo
2754          lsdistupdate = nldistupdate
2755
2756    END SELECT
2757
2758
2759 END SUBROUTINE set_salsa_runtime 
2760 
2761!------------------------------------------------------------------------------!
2762! Description:
2763! ------------
2764!> Calculates the absolute temperature (using hydrostatic pressure), saturation
2765!> vapour pressure and mixing ratio over water, relative humidity and air
2766!> density needed in the SALSA model.
2767!> NOTE, no saturation adjustment takes place -> the resulting water vapour
2768!> mixing ratio can be supersaturated, allowing the microphysical calculations
2769!> in SALSA.
2770!
2771!> Juha Tonttila, FMI, 2014 (original SALSAthrm)
2772!> Mona Kurppa, UHel, 2017 (adjustment for PALM and only aerosol processes)
2773!------------------------------------------------------------------------------!
2774 SUBROUTINE salsa_thrm_ij( i, j, p_ij, temp_ij, cw_ij, cs_ij, adn_ij )
2775 
2776    USE arrays_3d,                                                             &
2777        ONLY: p, pt, q, zu
2778       
2779    USE basic_constants_and_equations_mod,                                     &
2780        ONLY:  barometric_formula, exner_function, ideal_gas_law_rho, magnus
2781       
2782    USE control_parameters,                                                    &
2783        ONLY: pt_surface, surface_pressure
2784       
2785    IMPLICIT NONE
2786   
2787    INTEGER(iwp), INTENT(in) ::  i
2788    INTEGER(iwp), INTENT(in) ::  j
2789    REAL(wp), DIMENSION(:), INTENT(inout) ::  adn_ij
2790    REAL(wp), DIMENSION(:), INTENT(inout) ::  p_ij       
2791    REAL(wp), DIMENSION(:), INTENT(inout) ::  temp_ij
2792    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cw_ij
2793    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cs_ij
2794    REAL(wp), DIMENSION(nzb:nzt+1) ::  e_s !< saturation vapour pressure
2795                                           !< over water (Pa)
2796    REAL(wp) ::  t_surface !< absolute surface temperature (K)
2797!
2798!-- Pressure p_ijk (Pa) = hydrostatic pressure + perturbation pressure (p)
2799    t_surface = pt_surface * exner_function( surface_pressure * 100.0_wp )
2800    p_ij(:) = barometric_formula( zu, t_surface, surface_pressure * 100.0_wp ) &
2801              + p(:,j,i)
2802!             
2803!-- Absolute ambient temperature (K)
2804    temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) )       
2805!
2806!-- Air density
2807    adn_ij(:) = ideal_gas_law_rho( p_ij(:), temp_ij(:) )
2808!
2809!-- Water vapour concentration r_v (kg/m3)
2810    IF ( PRESENT( cw_ij ) )  THEN
2811       cw_ij(:) = ( q(:,j,i) / ( 1.0_wp - q(:,j,i) ) ) * adn_ij(:) 
2812    ENDIF
2813!
2814!-- Saturation mixing ratio r_s (kg/kg) from vapour pressure at temp (Pa)
2815    IF ( PRESENT( cs_ij ) )  THEN
2816       e_s(:) = 611.0_wp * EXP( alv_d_rv * ( 3.6609E-3_wp - 1.0_wp /           &
2817                temp_ij(:) ) )! magnus( temp_ij(:) )
2818       cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:) - e_s(:) ) ) * adn_ij(:) 
2819    ENDIF
2820
2821 END SUBROUTINE salsa_thrm_ij 
2822
2823!------------------------------------------------------------------------------!
2824! Description:
2825! ------------
2826!> Calculates ambient sizes of particles by equilibrating soluble fraction of
2827!> particles with water using the ZSR method (Stokes and Robinson, 1966).
2828!> Method:
2829!> Following chemical components are assumed water-soluble
2830!> - (ammonium) sulphate (100%)
2831!> - sea salt (100 %)
2832!> - organic carbon (epsoc * 100%)
2833!> Exact thermodynamic considerations neglected.
2834!> - If particles contain no sea salt, calculation according to sulphate
2835!>   properties
2836!> - If contain sea salt but no sulphate, calculation according to sea salt
2837!>   properties
2838!> - If contain both sulphate and sea salt -> the molar fraction of these
2839!>   compounds determines which one of them is used as the basis of calculation.
2840!> If sulphate and sea salt coexist in a particle, it is assumed that the Cl is
2841!> replaced by sulphate; thus only either sulphate + organics or sea salt +
2842!> organics is included in the calculation of soluble fraction.
2843!> Molality parameterizations taken from Table 1 of Tang: Thermodynamic and
2844!> optical properties of mixed-salt aerosols of atmospheric importance,
2845!> J. Geophys. Res., 102 (D2), 1883-1893 (1997)
2846!
2847!> Coded by:
2848!> Hannele Korhonen (FMI) 2005
2849!> Harri Kokkola (FMI) 2006
2850!> Matti Niskanen(FMI) 2012
2851!> Anton Laakso  (FMI) 2013
2852!> Modified for the new aerosol datatype, Juha Tonttila (FMI) 2014
2853!
2854!> fxm: should sea salt form a solid particle when prh is very low (even though
2855!> it could be mixed with e.g. sulphate)?
2856!> fxm: crashes if no sulphate or sea salt
2857!> fxm: do we really need to consider Kelvin effect for subrange 2
2858!------------------------------------------------------------------------------!     
2859 SUBROUTINE equilibration( prh, ptemp, paero, init )
2860     
2861    IMPLICIT NONE
2862!
2863!-- Input variables
2864    LOGICAL, INTENT(in) ::  init   !< TRUE: Initialization call
2865                                   !< FALSE: Normal runtime: update water
2866                                   !<        content only for 1a
2867    REAL(wp), INTENT(in) ::  prh   !< relative humidity [0-1]
2868    REAL(wp), INTENT(in) ::  ptemp !< temperature (K)
2869!
2870!-- Output variables
2871    TYPE(t_section), INTENT(inout) ::  paero(fn2b)     
2872!
2873!-- Local
2874    INTEGER(iwp) :: b      !< loop index
2875    INTEGER(iwp) :: counti  !< loop index
2876    REAL(wp) ::  zaw        !< water activity [0-1]       
2877    REAL(wp) ::  zbinmol(7) !< binary molality of each components (mol/kg)
2878    REAL(wp) ::  zcore      !< Volume of dry particle   
2879    REAL(wp) ::  zdold      !< Old diameter
2880    REAL(wp) ::  zdwet      !< Wet diameter or mean droplet diameter
2881    REAL(wp) ::  zke        !< Kelvin term in the Köhler equation
2882    REAL(wp) ::  zlwc       !< liquid water content [kg/m3-air]
2883    REAL(wp) ::  zrh        !< Relative humidity
2884    REAL(wp) ::  zvpart(7)  !< volume of chem. compounds in one particle
2885   
2886    zaw       = 0.0_wp
2887    zbinmol   = 0.0_wp
2888    zcore     = 0.0_wp
2889    zdold     = 0.0_wp
2890    zdwet     = 0.0_wp
2891    zlwc      = 0.0_wp
2892    zrh       = 0.0_wp
2893   
2894!               
2895!-- Relative humidity:
2896    zrh = prh
2897    zrh = MAX( zrh, 0.05_wp )
2898    zrh = MIN( zrh, 0.98_wp)   
2899!
2900!-- 1) Regime 1: sulphate and partly water-soluble OC. Done for every CALL
2901    DO  b = in1a, fn1a   ! size bin
2902         
2903       zbinmol = 0.0_wp
2904       zdold   = 1.0_wp 
2905       zke     = 1.02_wp
2906       
2907       IF ( paero(b)%numc > nclim )  THEN
2908!
2909!--       Volume in one particle
2910          zvpart = 0.0_wp
2911          zvpart(1:2) = paero(b)%volc(1:2) / paero(b)%numc
2912          zvpart(6:7) = paero(b)%volc(6:7) / paero(b)%numc
2913!               
2914!--       Total volume and wet diameter of one dry particle
2915          zcore = SUM( zvpart(1:2) )
2916          zdwet = paero(b)%dwet
2917         
2918          counti = 0
2919          DO  WHILE ( ABS( zdwet / zdold - 1.0_wp ) > 1.0E-2_wp ) 
2920         
2921             zdold = MAX( zdwet, 1.0E-20_wp )
2922             zaw = MAX( 1.0E-3_wp, zrh / zke ) ! To avoid underflow
2923!                   
2924!--          Binary molalities (mol/kg):
2925!--          Sulphate
2926             zbinmol(1) = 1.1065495E+2_wp - 3.6759197E+2_wp * zaw              &
2927                                          + 5.0462934E+2_wp * zaw**2.0_wp      &
2928                                          - 3.1543839E+2_wp * zaw**3.0_wp      &
2929                                          + 6.770824E+1_wp  * zaw**4.0_wp 
2930!--          Organic carbon                     
2931             zbinmol(2) = 1.0_wp / ( zaw * amh2o ) - 1.0_wp / amh2o 
2932!--          Nitric acid                             
2933             zbinmol(6) = 2.306844303E+1_wp - 3.563608869E+1_wp * zaw          &
2934                                            - 6.210577919E+1_wp * zaw**2.0_wp  &
2935                                            + 5.510176187E+2_wp * zaw**3.0_wp  &
2936                                            - 1.460055286E+3_wp * zaw**4.0_wp  &
2937                                            + 1.894467542E+3_wp * zaw**5.0_wp  &
2938                                            - 1.220611402E+3_wp * zaw**6.0_wp  &
2939                                            + 3.098597737E+2_wp * zaw**7.0_wp 
2940!
2941!--          Calculate the liquid water content (kg/m3-air) using ZSR (see e.g.
2942!--          Eq. 10.98 in Seinfeld and Pandis (2006))
2943             zlwc = ( paero(b)%volc(1) * ( arhoh2so4 / amh2so4 ) ) /           &
2944                    zbinmol(1) + epsoc * paero(b)%volc(2) * ( arhooc / amoc )  &
2945                    / zbinmol(2) + ( paero(b)%volc(6) * ( arhohno3/amhno3 ) )  &
2946                    / zbinmol(6)
2947!                           
2948!--          Particle wet diameter (m)
2949             zdwet = ( zlwc / paero(b)%numc / arhoh2o / api6 +                 &
2950                     ( SUM( zvpart(6:7) ) / api6 ) +      &
2951                       zcore / api6 )**( 1.0_wp / 3.0_wp )
2952!                             
2953!--          Kelvin effect (Eq. 10.85 in in Seinfeld and Pandis (2006)). Avoid
2954!--          overflow.
2955             zke = EXP( MIN( 50.0_wp,                                          &
2956                       4.0_wp * surfw0 * amvh2so4 / ( abo * ptemp *  zdwet ) ) )
2957             
2958             counti = counti + 1
2959             IF ( counti > 1000 )  THEN
2960                message_string = 'Subrange 1: no convergence!'
2961                CALL message( 'salsa_mod: equilibration', 'SA0042',            &
2962                              1, 2, 0, 6, 0 )
2963             ENDIF
2964          ENDDO
2965!               
2966!--       Instead of lwc, use the volume concentration of water from now on
2967!--       (easy to convert...)
2968          paero(b)%volc(8) = zlwc / arhoh2o
2969!               
2970!--       If this is initialization, update the core and wet diameter
2971          IF ( init )  THEN
2972             paero(b)%dwet = zdwet
2973             paero(b)%core = zcore
2974          ENDIF
2975         
2976       ELSE
2977!--       If initialization
2978!--       1.2) empty bins given bin average values 
2979          IF ( init )  THEN
2980             paero(b)%dwet = paero(b)%dmid
2981             paero(b)%core = api6 * paero(b)%dmid ** 3.0_wp
2982          ENDIF
2983         
2984       ENDIF
2985             
2986    ENDDO !< b
2987!
2988!-- 2) Regime 2a: sulphate, OC, BC and sea salt
2989!--    This is done only for initialization call, otherwise the water contents
2990!--    are computed via condensation
2991    IF ( init )  THEN
2992       DO  b = in2a, fn2b 
2993             
2994!--       Initialize
2995          zke     = 1.02_wp
2996          zbinmol = 0.0_wp
2997          zdold   = 1.0_wp
2998!               
2999!--       1) Particle properties calculated for non-empty bins
3000          IF ( paero(b)%numc > nclim )  THEN
3001!               
3002!--          Volume in one particle [fxm]
3003             zvpart = 0.0_wp
3004             zvpart(1:7) = paero(b)%volc(1:7) / paero(b)%numc
3005!
3006!--          Total volume and wet diameter of one dry particle [fxm]
3007             zcore = SUM( zvpart(1:5) )
3008             zdwet = paero(b)%dwet
3009
3010             counti = 0
3011             DO  WHILE ( ABS( zdwet / zdold - 1.0_wp ) > 1.0E-12_wp )
3012             
3013                zdold = MAX( zdwet, 1.0E-20_wp )
3014                zaw = zrh / zke
3015!                     
3016!--             Binary molalities (mol/kg):
3017!--             Sulphate
3018                zbinmol(1) = 1.1065495E+2_wp - 3.6759197E+2_wp * zaw           & 
3019                        + 5.0462934E+2_wp * zaw**2 - 3.1543839E+2_wp * zaw**3  &
3020                        + 6.770824E+1_wp  * zaw**4 
3021!--             Organic carbon                       
3022                zbinmol(2) = 1.0_wp / ( zaw * amh2o ) - 1.0_wp / amh2o 
3023!--             Nitric acid
3024                zbinmol(6) = 2.306844303E+1_wp - 3.563608869E+1_wp * zaw       &
3025                     - 6.210577919E+1_wp * zaw**2 + 5.510176187E+2_wp * zaw**3 &
3026                     - 1.460055286E+3_wp * zaw**4 + 1.894467542E+3_wp * zaw**5 &
3027                     - 1.220611402E+3_wp * zaw**6 + 3.098597737E+2_wp * zaw**7 
3028!--             Sea salt (natrium chloride)                                 
3029                zbinmol(5) = 5.875248E+1_wp - 1.8781997E+2_wp * zaw            &
3030                         + 2.7211377E+2_wp * zaw**2 - 1.8458287E+2_wp * zaw**3 &
3031                         + 4.153689E+1_wp  * zaw**4 
3032!                                 
3033!--             Calculate the liquid water content (kg/m3-air)
3034                zlwc = ( paero(b)%volc(1) * ( arhoh2so4 / amh2so4 ) ) /        &
3035                       zbinmol(1) + epsoc * ( paero(b)%volc(2) * ( arhooc /    &
3036                       amoc ) ) / zbinmol(2) + ( paero(b)%volc(6) * ( arhohno3 &
3037                       / amhno3 ) ) / zbinmol(6) + ( paero(b)%volc(5) *        &
3038                       ( arhoss / amss ) ) / zbinmol(5)
3039                       
3040!--             Particle wet radius (m)
3041                zdwet = ( zlwc / paero(b)%numc / arhoh2o / api6 +              &
3042                          ( SUM( zvpart(6:7) ) / api6 )  + &
3043                           zcore / api6 ) ** ( 1.0_wp / 3.0_wp )
3044!                               
3045!--             Kelvin effect (Eq. 10.85 in Seinfeld and Pandis (2006))
3046                zke = EXP( MIN( 50.0_wp,                                       &
3047                        4.0_wp * surfw0 * amvh2so4 / ( abo * zdwet * ptemp ) ) )
3048                         
3049                counti = counti + 1
3050                IF ( counti > 1000 )  THEN
3051                   message_string = 'Subrange 2: no convergence!'
3052                CALL message( 'salsa_mod: equilibration', 'SA0043',            &
3053                              1, 2, 0, 6, 0 )
3054                ENDIF
3055             ENDDO
3056!                   
3057!--          Liquid water content; instead of LWC use the volume concentration
3058             paero(b)%volc(8) = zlwc / arhoh2o
3059             paero(b)%dwet    = zdwet
3060             paero(b)%core    = zcore
3061             
3062          ELSE
3063!--          2.2) empty bins given bin average values
3064             paero(b)%dwet = paero(b)%dmid
3065             paero(b)%core = api6 * paero(b)%dmid ** 3.0_wp
3066          ENDIF
3067               
3068       ENDDO   ! b
3069    ENDIF
3070
3071 END SUBROUTINE equilibration 
3072 
3073!------------------------------------------------------------------------------!
3074!> Description:
3075!> ------------
3076!> Calculation of the settling velocity vc (m/s) per aerosol size bin and
3077!> deposition on plant canopy (lsdepo_vege).
3078!
3079!> Deposition is based on either the scheme presented in:
3080!> Zhang et al. (2001), Atmos. Environ. 35, 549-560 (includes collection due to
3081!> Brownian diffusion, impaction, interception and sedimentation)
3082!> OR
3083!> Petroff & Zhang (2010), Geosci. Model Dev. 3, 753-769 (includes also
3084!> collection due to turbulent impaction)
3085!
3086!> Equation numbers refer to equation in Jacobson (2005): Fundamentals of
3087!> Atmospheric Modeling, 2nd Edition.
3088!
3089!> Subroutine follows closely sedim_SALSA in UCLALES-SALSA written by Juha
3090!> Tonttila (KIT/FMI) and Zubair Maalick (UEF).
3091!> Rewritten to PALM by Mona Kurppa (UH), 2017.
3092!
3093!> Call for grid point i,j,k
3094!------------------------------------------------------------------------------!
3095
3096 SUBROUTINE deposition( paero, tk, adn, mag_u, lad, kvis, Sc, vc )
3097 
3098    USE plant_canopy_model_mod,                                                &
3099        ONLY: cdc
3100 
3101    IMPLICIT NONE
3102   
3103    REAL(wp), INTENT(in)    ::  adn    !< air density (kg/m3) 
3104    REAL(wp), INTENT(out)   ::  kvis   !< kinematic viscosity of air (m2/s)
3105    REAL(wp), INTENT(in) ::     lad    !< leaf area density (m2/m3)
3106    REAL(wp), INTENT(in)    ::  mag_u  !< wind velocity (m/s)
3107    REAL(wp), INTENT(out)   ::  Sc(:)  !< particle Schmidt number 
3108    REAL(wp), INTENT(in)    ::  tk     !< abs.temperature (K)   
3109    REAL(wp), INTENT(out)   ::  vc(:)  !< critical fall speed i.e. settling
3110                                       !< velocity of an aerosol particle (m/s)
3111    TYPE(t_section), INTENT(inout) ::  paero(fn2b)       
3112   
3113    INTEGER(iwp) ::  b      !< loop index
3114    REAL(wp) ::  avis       !< molecular viscocity of air (kg/(m*s))
3115    REAL(wp), PARAMETER ::  c_A = 1.249_wp !< Constants A, B and C for
3116    REAL(wp), PARAMETER ::  c_B = 0.42_wp  !< calculating  the Cunningham 
3117    REAL(wp), PARAMETER ::  c_C = 0.87_wp  !< slip-flow correction (Cc) 
3118                                           !< according to Jacobson (2005),
3119                                           !< Eq. 15.30
3120    REAL(wp) ::  Cc         !< Cunningham slip-flow correction factor     
3121    REAL(wp) ::  Kn         !< Knudsen number   
3122    REAL(wp) ::  lambda     !< molecular mean free path (m)
3123    REAL(wp) ::  mdiff      !< particle diffusivity coefficient   
3124    REAL(wp) ::  pdn        !< particle density (kg/m3)     
3125    REAL(wp) ::  ustar      !< friction velocity (m/s)   
3126    REAL(wp) ::  va         !< thermal speed of an air molecule (m/s)
3127    REAL(wp) ::  zdwet      !< wet diameter (m)                             
3128!
3129!-- Initialise
3130    Cc            = 0.0_wp
3131    Kn            = 0.0_wp
3132    mdiff         = 0.0_wp
3133    pdn           = 1500.0_wp    ! default value
3134    ustar         = 0.0_wp 
3135!
3136!-- Molecular viscosity of air (Eq. 4.54)
3137    avis = 1.8325E-5_wp * ( 416.16_wp / ( tk + 120.0_wp ) ) * ( tk /           &
3138           296.16_wp )**1.5_wp
3139!             
3140!-- Kinematic viscosity (Eq. 4.55)
3141    kvis =  avis / adn
3142!       
3143!-- Thermal velocity of an air molecule (Eq. 15.32)
3144    va = SQRT( 8.0_wp * abo * tk / ( pi * am_airmol ) ) 
3145!
3146!-- Mean free path (m) (Eq. 15.24)
3147    lambda = 2.0_wp * avis / ( adn * va )
3148   
3149    DO  b = 1, nbins
3150   
3151       IF ( paero(b)%numc < nclim )  CYCLE
3152       zdwet = paero(b)%dwet
3153!
3154!--    Knudsen number (Eq. 15.23)
3155       Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
3156!
3157!--    Cunningham slip-flow correction (Eq. 15.30)
3158       Cc = 1.0_wp + Kn * ( c_A + c_B * EXP( -c_C / Kn ) )
3159
3160!--    Particle diffusivity coefficient (Eq. 15.29)
3161       mdiff = ( abo * tk * Cc ) / ( 3.0_wp * pi * avis * zdwet )
3162!       
3163!--    Particle Schmidt number (Eq. 15.36)
3164       Sc(b) = kvis / mdiff       
3165!       
3166!--    Critical fall speed i.e. settling velocity  (Eq. 20.4)                 
3167       vc(b) = MIN( 1.0_wp, terminal_vel( 0.5_wp * zdwet, pdn, adn, avis, Cc) )
3168       
3169       IF ( lsdepo_vege  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
3170!       
3171!--       Friction velocity calculated following Prandtl (1925):
3172          ustar = SQRT( cdc ) * mag_u
3173          CALL depo_vege( paero, b, vc(b), mag_u, ustar, kvis, Sc(b), lad )
3174       ENDIF
3175    ENDDO
3176 
3177 END SUBROUTINE deposition 
3178 
3179!------------------------------------------------------------------------------!
3180! Description:
3181! ------------
3182!> Calculate change in number and volume concentrations due to deposition on
3183!> plant canopy.
3184!------------------------------------------------------------------------------!
3185 SUBROUTINE depo_vege( paero, b, vc, mag_u, ustar, kvis_a, Sc, lad )
3186 
3187    IMPLICIT NONE
3188   
3189    INTEGER(iwp), INTENT(in) ::  b  !< loop index
3190    REAL(wp), INTENT(in) ::  kvis_a !< kinematic viscosity of air (m2/s)
3191    REAL(wp), INTENT(in) ::  lad    !< leaf area density (m2/m3)
3192    REAL(wp), INTENT(in) ::  mag_u  !< wind velocity (m/s)   
3193    REAL(wp), INTENT(in) ::  Sc     !< particle Schmidt number
3194    REAL(wp), INTENT(in) ::  ustar  !< friction velocity (m/s)                                   
3195    REAL(wp), INTENT(in) ::  vc     !< terminal velocity (m/s) 
3196    TYPE(t_section), INTENT(inout) ::  paero(fn2b) 
3197   
3198    INTEGER(iwp) ::  c      !< loop index
3199
3200    REAL(wp) ::  alpha       !< parameter, Table 3 in Zhang et al. (2001) 
3201    REAL(wp) ::  depo        !< deposition efficiency
3202    REAL(wp) ::  C_Br        !< coefficient for Brownian diffusion
3203    REAL(wp) ::  C_IM        !< coefficient for inertial impaction
3204    REAL(wp) ::  C_IN        !< coefficient for interception
3205    REAL(wp) ::  C_IT        !< coefficient for turbulent impaction   
3206    REAL(wp) ::  gamma       !< parameter, Table 3 in Zhang et al. (2001)   
3207    REAL(wp) ::  par_A       !< parameter A for the characteristic radius of
3208                             !< collectors, Table 3 in Zhang et al. (2001)   
3209    REAL(wp) ::  rt          !< the overall quasi-laminar resistance for
3210                             !< particles
3211    REAL(wp) ::  St          !< Stokes number for smooth surfaces or bluff
3212                             !< surface elements                                 
3213    REAL(wp) ::  tau_plus    !< dimensionless particle relaxation time   
3214    REAL(wp) ::  v_bd        !< deposition velocity due to Brownian diffusion
3215    REAL(wp) ::  v_im        !< deposition velocity due to impaction
3216    REAL(wp) ::  v_in        !< deposition velocity due to interception
3217    REAL(wp) ::  v_it        !< deposition velocity due to turbulent impaction                               
3218!
3219!-- Initialise
3220    depo     = 0.0_wp 
3221    rt       = 0.0_wp
3222    St       = 0.0_wp
3223    tau_plus = 0.0_wp
3224    v_bd     = 0.0_wp     
3225    v_im     = 0.0_wp       
3226    v_in     = 0.0_wp       
3227    v_it     = 0.0_wp         
3228       
3229    IF ( depo_vege_type == 'zhang2001' )  THEN
3230!       
3231!--    Parameters for the land use category 'deciduous broadleaf trees'(Table 3)     
3232       par_A = 5.0E-3_wp
3233       alpha = 0.8_wp
3234       gamma = 0.56_wp 
3235!       
3236!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24) 
3237       St = vc * ustar / ( g * par_A )         
3238!         
3239!--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)       
3240       rt = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -St**0.5_wp ) *    &
3241                         ( Sc**( -gamma ) + ( St / ( alpha + St ) )**2.0_wp +  &
3242                           0.5_wp * ( paero(b)%dwet / par_A )**2.0_wp ) ) )
3243       depo = ( rt + vc ) * lad
3244       paero(b)%numc = paero(b)%numc - depo * paero(b)%numc * dt_salsa
3245       DO  c = 1, maxspec+1
3246          paero(b)%volc(c) = paero(b)%volc(c) - depo * paero(b)%volc(c) *      &
3247                             dt_salsa
3248       ENDDO
3249       
3250    ELSEIF ( depo_vege_type == 'petroff2010' )  THEN
3251!
3252!--    vd = v_BD + v_IN + v_IM + v_IT + vc
3253!--    Deposition efficiencies from Table 1. Constants from Table 2.
3254       C_Br  = 1.262_wp
3255       C_IM  = 0.130_wp
3256       C_IN  = 0.216_wp
3257       C_IT  = 0.056_wp
3258       par_A = 0.03_wp   ! Here: leaf width (m)     
3259!       
3260!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24) 
3261       St = vc * ustar / ( g * par_A )         
3262!
3263!--    Non-dimensional relexation time of the particle on top of canopy
3264       tau_plus = vc * ustar**2.0_wp / ( kvis_a * g ) 
3265!
3266!--    Brownian diffusion
3267       v_bd = mag_u * C_Br * Sc**( -2.0_wp / 3.0_wp ) *                        &
3268              ( mag_u * par_A / kvis_a )**( -0.5_wp )
3269!
3270!--    Interception
3271       v_in = mag_u * C_IN * paero(b)%dwet / par_A * ( 2.0_wp + LOG( 2.0_wp *  &
3272              par_A / paero(b)%dwet ) )                     
3273!
3274!--    Impaction: Petroff (2009) Eq. 18
3275       v_im = mag_u * C_IM * ( St / ( St + 0.47_wp ) )**2.0_wp
3276       
3277       IF ( tau_plus < 20.0_wp )  THEN
3278          v_it = 2.5E-3_wp * C_IT * tau_plus**2.0_wp
3279       ELSE
3280          v_it = C_IT
3281       ENDIF
3282       depo = ( v_bd + v_in + v_im + v_it + vc ) * lad     
3283       paero(b)%numc = paero(b)%numc - depo * paero(b)%numc * dt_salsa     
3284       DO  c = 1, maxspec+1
3285          paero(b)%volc(c) = paero(b)%volc(c) - depo * paero(b)%volc(c) *      &
3286                             dt_salsa
3287       ENDDO
3288    ENDIF 
3289 
3290 END SUBROUTINE depo_vege
3291 
3292!------------------------------------------------------------------------------!
3293! Description:
3294! ------------ 
3295!> Calculate deposition on horizontal and vertical surfaces. Implement as
3296!> surface flux.
3297!------------------------------------------------------------------------------!
3298
3299 SUBROUTINE depo_topo( i, j, surf, vc, Sc, kvis, mag_u, norm )
3300 
3301    USE surface_mod,                                                           &
3302        ONLY:  surf_type
3303 
3304    IMPLICIT NONE
3305   
3306    INTEGER(iwp), INTENT(in) ::  i     !< loop index
3307    INTEGER(iwp), INTENT(in) ::  j     !< loop index
3308    REAL(wp), INTENT(in) ::  kvis(:)   !< kinematic viscosity of air (m2/s)
3309    REAL(wp), INTENT(in) ::  mag_u(:)  !< wind velocity (m/s)                                                 
3310    REAL(wp), INTENT(in) ::  norm(:)   !< normalisation (usually air density)
3311    REAL(wp), INTENT(in) ::  Sc(:,:)  !< particle Schmidt number
3312    REAL(wp), INTENT(in) ::  vc(:,:)  !< terminal velocity (m/s)   
3313    TYPE(surf_type), INTENT(inout) :: surf  !< respective surface type
3314    INTEGER(iwp) ::  b      !< loop index
3315    INTEGER(iwp) ::  c      !< loop index
3316    INTEGER(iwp) ::  k      !< loop index
3317    INTEGER(iwp) ::  m      !< loop index
3318    INTEGER(iwp) ::  surf_e !< End index of surface elements at (j,i)-gridpoint
3319    INTEGER(iwp) ::  surf_s !< Start index of surface elements at (j,i)-gridpoint
3320    REAL(wp) ::  alpha      !< parameter, Table 3 in Zhang et al. (2001)
3321    REAL(wp) ::  C_Br       !< coefficient for Brownian diffusion
3322    REAL(wp) ::  C_IM       !< coefficient for inertial impaction
3323    REAL(wp) ::  C_IN       !< coefficient for interception
3324    REAL(wp) ::  C_IT       !< coefficient for turbulent impaction
3325    REAL(wp) ::  depo       !< deposition efficiency
3326    REAL(wp) ::  gamma      !< parameter, Table 3 in Zhang et al. (2001)
3327    REAL(wp) ::  par_A      !< parameter A for the characteristic radius of
3328                            !< collectors, Table 3 in Zhang et al. (2001)
3329    REAL(wp) ::  rt         !< the overall quasi-laminar resistance for
3330                            !< particles
3331    REAL(wp) ::  St         !< Stokes number for bluff surface elements 
3332    REAL(wp) ::  tau_plus   !< dimensionless particle relaxation time   
3333    REAL(wp) ::  v_bd       !< deposition velocity due to Brownian diffusion
3334    REAL(wp) ::  v_im       !< deposition velocity due to impaction
3335    REAL(wp) ::  v_in       !< deposition velocity due to interception
3336    REAL(wp) ::  v_it       !< deposition velocity due to turbulent impaction 
3337!
3338!-- Initialise
3339    rt       = 0.0_wp
3340    St       = 0.0_wp
3341    tau_plus = 0.0_wp
3342    v_bd     = 0.0_wp     
3343    v_im     = 0.0_wp       
3344    v_in     = 0.0_wp       
3345    v_it     = 0.0_wp                                 
3346    surf_s   = surf%start_index(j,i)
3347    surf_e   = surf%end_index(j,i) 
3348   
3349    DO  m = surf_s, surf_e 
3350       k = surf%k(m)       
3351       DO  b = 1, nbins
3352          IF ( aerosol_number(b)%conc(k,j,i) <= nclim  .OR.                    &
3353               Sc(k+1,b) < 1.0_wp )  CYCLE   
3354                   
3355          IF ( depo_topo_type == 'zhang2001' )  THEN
3356!       
3357!--          Parameters for the land use category 'urban' in Table 3
3358             alpha = 1.5_wp
3359             gamma = 0.56_wp 
3360             par_A = 10.0E-3_wp
3361!       
3362!--          Stokes number for smooth surfaces or surfaces with bluff roughness
3363!--          elements (Seinfeld and Pandis, 2nd edition (2006): Eq. 19.23)       
3364             St = MAX( 0.01_wp, vc(k+1,b) * surf%us(m) ** 2.0_wp /             &
3365                       ( g * kvis(k+1)  ) ) 
3366!         
3367!--          The overall quasi-laminar resistance for particles (Eq. 5)       
3368             rt = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * surf%us(m) * (            &
3369                       Sc(k+1,b)**( -gamma ) + ( St / ( alpha + St ) )**2.0_wp &
3370                        + 0.5_wp * ( Ra_dry(k,j,i,b) / par_A )**2.0_wp ) *     &
3371                       EXP( -St**0.5_wp ) ) ) 
3372             depo = vc(k+1,b) + rt
3373             
3374          ELSEIF ( depo_topo_type == 'petroff2010' )  THEN 
3375!
3376!--          vd = v_BD + v_IN + v_IM + v_IT + vc
3377!--          Deposition efficiencies from Table 1. Constants from Table 2.
3378             C_Br  = 1.262_wp
3379             C_IM  = 0.130_wp
3380             C_IN  = 0.216_wp
3381             C_IT  = 0.056_wp
3382             par_A = 0.03_wp   ! Here: leaf width (m) 
3383!       
3384!--          Stokes number for smooth surfaces or surfaces with bluff roughness
3385!--          elements (Seinfeld and Pandis, 2nd edition (2006): Eq. 19.23)       
3386             St = MAX( 0.01_wp, vc(k+1,b) * surf%us(m) ** 2.0_wp /             &
3387                       ( g *  kvis(k+1) ) )             
3388!
3389!--          Non-dimensional relexation time of the particle on top of canopy
3390             tau_plus = vc(k+1,b) * surf%us(m)**2.0_wp / ( kvis(k+1) * g ) 
3391!
3392!--          Brownian diffusion
3393             v_bd = mag_u(k+1) * C_Br * Sc(k+1,b)**( -2.0_wp / 3.0_wp ) *      &
3394                    ( mag_u(k+1) * par_A / kvis(k+1) )**( -0.5_wp )
3395!
3396!--          Interception
3397             v_in = mag_u(k+1) * C_IN * Ra_dry(k,j,i,b)/ par_A * ( 2.0_wp +    &
3398                    LOG( 2.0_wp * par_A / Ra_dry(k,j,i,b) ) )                     
3399!
3400!--          Impaction: Petroff (2009) Eq. 18
3401             v_im = mag_u(k+1) * C_IM * ( St / ( St + 0.47_wp ) )**2.0_wp
3402             
3403             IF ( tau_plus < 20.0_wp )  THEN
3404                v_it = 2.5E-3_wp * C_IT * tau_plus**2.0_wp
3405             ELSE
3406                v_it = C_IT
3407             ENDIF
3408             depo =  v_bd + v_in + v_im + v_it + vc(k+1,b)       
3409         
3410          ENDIF
3411          IF ( lod_aero == 3  .OR.  salsa_source_mode ==  'no_source' )  THEN
3412             surf%answs(m,b) = -depo * norm(k) * aerosol_number(b)%conc(k,j,i) 
3413             DO  c = 1, ncc_tot   
3414                surf%amsws(m,(c-1)*nbins+b) = -depo *  norm(k) *               &
3415                                         aerosol_mass((c-1)*nbins+b)%conc(k,j,i)
3416             ENDDO    ! c
3417          ELSE
3418             surf%answs(m,b) = SUM( aerosol_number(b)%source(:,j,i) ) -        &
3419                               MAX( 0.0_wp, depo * norm(k) *                   &
3420                               aerosol_number(b)%conc(k,j,i) )
3421             DO  c = 1, ncc_tot   
3422                surf%amsws(m,(c-1)*nbins+b) = SUM(                             &
3423                               aerosol_mass((c-1)*nbins+b)%source(:,j,i) ) -   &
3424                               MAX(  0.0_wp, depo *  norm(k) *                 &
3425                               aerosol_mass((c-1)*nbins+b)%conc(k,j,i) )
3426             ENDDO 
3427          ENDIF
3428       ENDDO    ! b
3429    ENDDO    ! m     
3430     
3431 END SUBROUTINE depo_topo
3432 
3433!------------------------------------------------------------------------------!
3434! Description:
3435! ------------
3436! Function for calculating terminal velocities for different particles sizes.
3437!------------------------------------------------------------------------------!
3438 REAL(wp) FUNCTION terminal_vel( radius, rhop, rhoa, visc, beta )
3439 
3440    IMPLICIT NONE
3441   
3442    REAL(wp), INTENT(in) ::  beta    !< Cunningham correction factor
3443    REAL(wp), INTENT(in) ::  radius  !< particle radius (m)
3444    REAL(wp), INTENT(in) ::  rhop    !< particle density (kg/m3)
3445    REAL(wp), INTENT(in) ::  rhoa    !< air density (kg/m3)
3446    REAL(wp), INTENT(in) ::  visc    !< molecular viscosity of air (kg/(m*s))
3447   
3448!
3449!-- Stokes law with Cunningham slip correction factor
3450    terminal_vel = ( 4.0_wp * radius**2.0_wp ) * ( rhop - rhoa ) * g * beta /  &
3451                   ( 18.0_wp * visc ) ! (m/s)
3452       
3453 END FUNCTION terminal_vel
3454 
3455!------------------------------------------------------------------------------!
3456! Description:
3457! ------------
3458!> Calculates particle loss and change in size distribution due to (Brownian)
3459!> coagulation. Only for particles with dwet < 30 micrometres.
3460!
3461!> Method:
3462!> Semi-implicit, non-iterative method: (Jacobson, 1994)
3463!> Volume concentrations of the smaller colliding particles added to the bin of
3464!> the larger colliding particles. Start from first bin and use the updated
3465!> number and volume for calculation of following bins. NB! Our bin numbering
3466!> does not follow particle size in subrange 2.
3467!
3468!> Schematic for bin numbers in different subranges:
3469!>             1                            2
3470!>    +-------------------------------------------+
3471!>  a | 1 | 2 | 3 || 4 | 5 | 6 | 7 |  8 |  9 | 10||
3472!>  b |           ||11 |12 |13 |14 | 15 | 16 | 17||
3473!>    +-------------------------------------------+
3474!
3475!> Exact coagulation coefficients for each pressure level are scaled according
3476!> to current particle wet size (linear scaling).
3477!> Bins are organized in terms of the dry size of the condensation nucleus,
3478!> while coagulation kernell is calculated with the actual hydrometeor
3479!> size.
3480!
3481!> Called from salsa_driver
3482!> fxm: Process selection should be made smarter - now just lots of IFs inside
3483!>      loops
3484!
3485!> Coded by:
3486!> Hannele Korhonen (FMI) 2005
3487!> Harri Kokkola (FMI) 2006
3488!> Tommi Bergman (FMI) 2012
3489!> Matti Niskanen(FMI) 2012
3490!> Anton Laakso  (FMI) 2013
3491!> Juha Tonttila (FMI) 2014
3492!------------------------------------------------------------------------------!
3493 SUBROUTINE coagulation( paero, ptstep, ptemp, ppres )
3494               
3495    IMPLICIT NONE
3496   
3497!-- Input and output variables
3498    TYPE(t_section), INTENT(inout) ::  paero(fn2b) !< Aerosol properties
3499    REAL(wp), INTENT(in) ::  ppres  !< ambient pressure (Pa)
3500    REAL(wp), INTENT(in) ::  ptemp  !< ambient temperature (K)
3501    REAL(wp), INTENT(in) ::  ptstep !< time step (s)
3502!-- Local variables
3503    INTEGER(iwp) ::  index_2a !< corresponding bin in subrange 2a
3504    INTEGER(iwp) ::  index_2b !< corresponding bin in subrange 2b
3505    INTEGER(iwp) ::  b !< loop index
3506    INTEGER(iwp) ::  ll !< loop index
3507    INTEGER(iwp) ::  mm !< loop index
3508    INTEGER(iwp) ::  nn !< loop index
3509    REAL(wp) ::  pressi !< pressure
3510    REAL(wp) ::  temppi !< temperature
3511    REAL(wp) ::  zcc(fn2b,fn2b)   !< updated coagulation coefficients (m3/s) 
3512    REAL(wp) ::  zdpart_mm        !< diameter of particle (m)
3513    REAL(wp) ::  zdpart_nn        !< diameter of particle (m)   
3514    REAL(wp) ::  zminusterm       !< coagulation loss in a bin (1/s)
3515    REAL(wp) ::  zplusterm(8)     !< coagulation gain in a bin (fxm/s)
3516                                  !< (for each chemical compound)
3517    REAL(wp) ::  zmpart(fn2b)     !< approximate mass of particles (kg)
3518   
3519    zcc       = 0.0_wp
3520    zmpart    = 0.0_wp
3521    zdpart_mm = 0.0_wp
3522    zdpart_nn = 0.0_wp
3523!
3524!-- 1) Coagulation to coarse mode calculated in a simplified way:
3525!--    CoagSink ~ Dp in continuum subrange, thus we calculate 'effective'
3526!--    number concentration of coarse particles
3527
3528!-- 2) Updating coagulation coefficients
3529!   
3530!-- Aerosol mass (kg). Density of 1500 kg/m3 assumed
3531    zmpart(1:fn2b) = api6 * ( MIN( paero(1:fn2b)%dwet, 30.0E-6_wp )**3.0_wp  ) &
3532                     * 1500.0_wp 
3533    temppi = ptemp
3534    pressi = ppres
3535    zcc    = 0.0_wp
3536!
3537!-- Aero-aero coagulation
3538    DO  mm = 1, fn2b   ! smaller colliding particle
3539       IF ( paero(mm)%numc < nclim )  CYCLE
3540       DO  nn = mm, fn2b   ! larger colliding particle
3541          IF ( paero(nn)%numc < nclim )  CYCLE
3542         
3543          zdpart_mm = MIN( paero(mm)%dwet, 30.0E-6_wp )     ! Limit to 30 um
3544          zdpart_nn = MIN( paero(nn)%dwet, 30.0E-6_wp )     ! Limit to 30 um
3545!             
3546!--       Coagulation coefficient of particles (m3/s)
3547          zcc(mm,nn) = coagc( zdpart_mm, zdpart_nn, zmpart(mm), zmpart(nn),    &
3548                              temppi, pressi )
3549          zcc(nn,mm) = zcc(mm,nn)
3550       ENDDO
3551    ENDDO
3552       
3553!   
3554!-- 3) New particle and volume concentrations after coagulation:
3555!--    Calculated according to Jacobson (2005) eq. 15.9
3556!
3557!-- Aerosols in subrange 1a:
3558    DO  b = in1a, fn1a
3559       IF ( paero(b)%numc < nclim )  CYCLE
3560       zminusterm   = 0.0_wp
3561       zplusterm(:) = 0.0_wp
3562!       
3563!--    Particles lost by coagulation with larger aerosols
3564       DO  ll = b+1, fn2b
3565          zminusterm = zminusterm + zcc(b,ll) * paero(ll)%numc
3566       ENDDO
3567!       
3568!--    Coagulation gain in a bin: change in volume conc. (cm3/cm3):
3569       DO ll = in1a, b-1
3570          zplusterm(1:2) = zplusterm(1:2) + zcc(ll,b) * paero(ll)%volc(1:2)
3571          zplusterm(6:7) = zplusterm(6:7) + zcc(ll,b) * paero(ll)%volc(6:7)
3572          zplusterm(8)   = zplusterm(8)   + zcc(ll,b) * paero(ll)%volc(8)
3573       ENDDO
3574!       
3575!--    Volume and number concentrations after coagulation update [fxm]
3576       paero(b)%volc(1:2) = ( paero(b)%volc(1:2) + ptstep * zplusterm(1:2) * &
3577                             paero(b)%numc ) / ( 1.0_wp + ptstep * zminusterm )
3578       paero(b)%volc(6:7) = ( paero(b)%volc(6:7) + ptstep * zplusterm(6:7) * &
3579                             paero(b)%numc ) / ( 1.0_wp + ptstep * zminusterm )
3580       paero(b)%volc(8)   = ( paero(b)%volc(8)   + ptstep * zplusterm(8) *   &
3581                             paero(b)%numc ) / ( 1.0_wp + ptstep * zminusterm )
3582       paero(b)%numc = paero(b)%numc / ( 1.0_wp + ptstep * zminusterm  +     &
3583                        0.5_wp * ptstep * zcc(b,b) * paero(b)%numc )               
3584    ENDDO
3585!             
3586!-- Aerosols in subrange 2a:
3587    DO  b = in2a, fn2a
3588       IF ( paero(b)%numc < nclim )  CYCLE
3589       zminusterm   = 0.0_wp
3590       zplusterm(:) = 0.0_wp
3591!       
3592!--    Find corresponding size bin in subrange 2b
3593       index_2b = b - in2a + in2b
3594!       
3595!--    Particles lost by larger particles in 2a
3596       DO  ll = b+1, fn2a
3597          zminusterm = zminusterm + zcc(b,ll) * paero(ll)%numc 
3598       ENDDO
3599!       
3600!--    Particles lost by larger particles in 2b
3601       IF ( .NOT. no_insoluble )  THEN
3602          DO  ll = index_2b+1, fn2b
3603             zminusterm = zminusterm + zcc(b,ll) * paero(ll)%numc
3604          ENDDO
3605       ENDIF
3606!       
3607!--    Particle volume gained from smaller particles in subranges 1, 2a and 2b
3608       DO  ll = in1a, b-1
3609          zplusterm(1:2) = zplusterm(1:2) + zcc(ll,b) * paero(ll)%volc(1:2)
3610          zplusterm(6:7) = zplusterm(6:7) + zcc(ll,b) * paero(ll)%volc(6:7)
3611          zplusterm(8)   = zplusterm(8)   + zcc(ll,b) * paero(ll)%volc(8)
3612       ENDDO 
3613!       
3614!--    Particle volume gained from smaller particles in 2a
3615!--    (Note, for components not included in the previous loop!)
3616       DO  ll = in2a, b-1
3617          zplusterm(3:5) = zplusterm(3:5) + zcc(ll,b)*paero(ll)%volc(3:5)             
3618       ENDDO
3619       
3620!       
3621!--    Particle volume gained from smaller (and equal) particles in 2b
3622       IF ( .NOT. no_insoluble )  THEN
3623          DO  ll = in2b, index_2b
3624             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,b) * paero(ll)%volc(1:8)
3625          ENDDO
3626       ENDIF
3627!       
3628!--    Volume and number concentrations after coagulation update [fxm]
3629       paero(b)%volc(1:8) = ( paero(b)%volc(1:8) + ptstep * zplusterm(1:8) * &
3630                             paero(b)%numc ) / ( 1.0_wp + ptstep * zminusterm )
3631       paero(b)%numc = paero(b)%numc / ( 1.0_wp + ptstep * zminusterm +      &
3632                        0.5_wp * ptstep * zcc(b,b) * paero(b)%numc )
3633    ENDDO
3634!             
3635!-- Aerosols in subrange 2b:
3636    IF ( .NOT. no_insoluble )  THEN
3637       DO  b = in2b, fn2b
3638          IF ( paero(b)%numc < nclim )  CYCLE
3639          zminusterm   = 0.0_wp
3640          zplusterm(:) = 0.0_wp
3641!       
3642!--       Find corresponding size bin in subsubrange 2a
3643          index_2a = b - in2b + in2a
3644!       
3645!--       Particles lost to larger particles in subranges 2b
3646          DO  ll = b+1, fn2b
3647             zminusterm = zminusterm + zcc(b,ll) * paero(ll)%numc
3648          ENDDO
3649!       
3650!--       Particles lost to larger and equal particles in 2a
3651          DO  ll = index_2a, fn2a
3652             zminusterm = zminusterm + zcc(b,ll) * paero(ll)%numc
3653          ENDDO
3654!       
3655!--       Particle volume gained from smaller particles in subranges 1 & 2a
3656          DO  ll = in1a, index_2a-1
3657             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,b) * paero(ll)%volc(1:8)
3658          ENDDO
3659!       
3660!--       Particle volume gained from smaller particles in 2b
3661          DO  ll = in2b, b-1
3662             zplusterm(1:8) = zplusterm(1:8) + zcc(ll,b) * paero(ll)%volc(1:8)
3663          ENDDO
3664!       
3665!--       Volume and number concentrations after coagulation update [fxm]
3666          paero(b)%volc(1:8) = ( paero(b)%volc(1:8) + ptstep * zplusterm(1:8)&
3667                           * paero(b)%numc ) / ( 1.0_wp + ptstep * zminusterm )
3668          paero(b)%numc = paero(b)%numc / ( 1.0_wp + ptstep * zminusterm  +  &
3669                           0.5_wp * ptstep * zcc(b,b) * paero(b)%numc )
3670       ENDDO
3671    ENDIF
3672
3673 END SUBROUTINE coagulation
3674
3675!------------------------------------------------------------------------------!
3676! Description:
3677! ------------
3678!> Calculation of coagulation coefficients. Extended version of the function
3679!> originally found in mo_salsa_init.
3680!
3681!> J. Tonttila, FMI, 05/2014
3682!------------------------------------------------------------------------------!
3683 REAL(wp) FUNCTION coagc( diam1, diam2, mass1, mass2, temp, pres )
3684 
3685    IMPLICIT NONE
3686!       
3687!-- Input and output variables
3688    REAL(wp), INTENT(in) ::  diam1 !< diameter of colliding particle 1 (m)
3689    REAL(wp), INTENT(in) ::  diam2 !< diameter of colliding particle 2 (m)
3690    REAL(wp), INTENT(in) ::  mass1 !< mass of colliding particle 1 (kg)
3691    REAL(wp), INTENT(in) ::  mass2 !< mass of colliding particle 2 (kg)
3692    REAL(wp), INTENT(in) ::  pres  !< ambient pressure (Pa?) [fxm]
3693    REAL(wp), INTENT(in) ::  temp  !< ambient temperature (K)       
3694!
3695!-- Local variables
3696    REAL(wp) ::  fmdist !< distance of flux matching (m)   
3697    REAL(wp) ::  knud_p !< particle Knudsen number
3698    REAL(wp) ::  mdiam  !< mean diameter of colliding particles (m) 
3699    REAL(wp) ::  mfp    !< mean free path of air molecules (m)   
3700    REAL(wp) ::  visc   !< viscosity of air (kg/(m s))                   
3701    REAL(wp), DIMENSION (2) ::  beta   !< Cunningham correction factor
3702    REAL(wp), DIMENSION (2) ::  dfpart !< particle diffusion coefficient
3703                                       !< (m2/s)       
3704    REAL(wp), DIMENSION (2) ::  diam   !< diameters of particles (m)
3705    REAL(wp), DIMENSION (2) ::  flux   !< flux in continuum and free molec.
3706                                       !< regime (m/s)       
3707    REAL(wp), DIMENSION (2) ::  knud   !< particle Knudsen number       
3708    REAL(wp), DIMENSION (2) ::  mpart  !< masses of particles (kg)
3709    REAL(wp), DIMENSION (2) ::  mtvel  !< particle mean thermal velocity (m/s)
3710    REAL(wp), DIMENSION (2) ::  omega  !< particle mean free path             
3711    REAL(wp), DIMENSION (2) ::  tva    !< temporary variable (m)       
3712!
3713!-- Initialisation
3714    coagc   = 0.0_wp
3715!
3716!-- 1) Initializing particle and ambient air variables
3717    diam  = (/ diam1, diam2 /) !< particle diameters (m)
3718    mpart = (/ mass1, mass2 /) !< particle masses (kg)
3719!-- Viscosity of air (kg/(m s))       
3720    visc = ( 7.44523E-3_wp * temp ** 1.5_wp ) /                                &
3721           ( 5093.0_wp * ( temp + 110.4_wp ) ) 
3722!-- Mean free path of air (m)           
3723    mfp = ( 1.656E-10_wp * temp + 1.828E-8_wp ) * ( p_0 + 1325.0_wp ) / pres
3724!
3725!-- 2) Slip correction factor for small particles
3726    knud = 2.0_wp * EXP( LOG(mfp) - LOG(diam) )! Knudsen number for air (15.23)
3727!-- Cunningham correction factor (Allen and Raabe, Aerosol Sci. Tech. 4, 269)       
3728    beta = 1.0_wp + knud * ( 1.142_wp + 0.558_wp * EXP( -0.999_wp / knud ) )
3729!
3730!-- 3) Particle properties
3731!-- Diffusion coefficient (m2/s) (Jacobson (2005) eq. 15.29)
3732    dfpart = beta * abo * temp / ( 3.0_wp * pi * visc * diam ) 
3733!-- Mean thermal velocity (m/s) (Jacobson (2005) eq. 15.32)
3734    mtvel = SQRT( ( 8.0_wp * abo * temp ) / ( pi * mpart ) )
3735!-- Particle mean free path (m) (Jacobson (2005) eq. 15.34 )
3736    omega = 8.0_wp * dfpart / ( pi * mtvel ) 
3737!-- Mean diameter (m)
3738    mdiam = 0.5_wp * ( diam(1) + diam(2) )
3739!
3740!-- 4) Calculation of fluxes (Brownian collision kernels) and flux matching
3741!-- following Jacobson (2005):
3742!-- Flux in continuum regime (m3/s) (eq. 15.28)
3743    flux(1) = 4.0_wp * pi * mdiam * ( dfpart(1) + dfpart(2) )
3744!-- Flux in free molec. regime (m3/s) (eq. 15.31)
3745    flux(2) = pi * SQRT( ( mtvel(1)**2.0_wp ) + ( mtvel(2)**2.0_wp ) ) *      &
3746              ( mdiam**2.0_wp )
3747!-- temporary variables (m) to calculate flux matching distance (m)
3748    tva(1) = ( ( mdiam + omega(1) )**3.0_wp - ( mdiam**2.0_wp +                &
3749               omega(1)**2.0_wp ) * SQRT( ( mdiam**2.0_wp + omega(1)**2.0_wp ) &
3750               ) ) / ( 3.0_wp * mdiam * omega(1) ) - mdiam
3751    tva(2) = ( ( mdiam + omega(2) )**3.0_wp - ( mdiam**2.0_wp +                &
3752               omega(2)**2.0_wp ) * SQRT( ( mdiam**2 + omega(2)**2 ) ) ) /     &
3753             ( 3.0_wp * mdiam * omega(2) ) - mdiam
3754!-- Flux matching distance (m) i.e. the mean distance from the centre of a
3755!-- sphere reached by particles leaving sphere's surface and travelling a
3756!-- distance of particle mean free path mfp (eq. 15 34)                 
3757    fmdist = SQRT( tva(1)**2 + tva(2)**2.0_wp) 
3758!
3759!-- 5) Coagulation coefficient (m3/s) (eq. 15.33). Here assumed
3760!-- coalescence efficiency 1!!
3761    coagc = flux(1) / ( mdiam / ( mdiam + fmdist) + flux(1) / flux(2) ) 
3762!-- coagulation coefficient = coalescence efficiency * collision kernel
3763!
3764!-- Corrected collision kernel following Karl et al., 2016 (ACP):
3765!-- Inclusion of van der Waals and viscous forces
3766    IF ( van_der_waals_coagc )  THEN
3767       knud_p = SQRT( omega(1)**2 + omega(2)**2 ) / mdiam   
3768       IF ( knud_p >= 0.1_wp  .AND.  knud_p <= 10.0_wp )  THEN
3769          coagc = coagc * ( 2.0_wp + 0.4_wp * LOG( knud_p ) )
3770       ELSE
3771          coagc = coagc * 3.0_wp
3772       ENDIF
3773    ENDIF
3774   
3775 END FUNCTION coagc
3776 
3777!------------------------------------------------------------------------------!   
3778! Description:
3779! ------------
3780!> Calculates the change in particle volume and gas phase
3781!> concentrations due to nucleation, condensation and dissolutional growth.
3782!
3783!> Sulphuric acid and organic vapour: only condensation and no evaporation.
3784!
3785!> New gas and aerosol phase concentrations calculated according to Jacobson
3786!> (1997): Numerical techniques to solve condensational and dissolutional growth
3787!> equations when growth is coupled to reversible reactions, Aerosol Sci. Tech.,
3788!> 27, pp 491-498.
3789!
3790!> Following parameterization has been used:
3791!> Molecular diffusion coefficient of condensing vapour (m2/s)
3792!> (Reid et al. (1987): Properties of gases and liquids, McGraw-Hill, New York.)
3793!> D = {1.d-7*sqrt(1/M_air + 1/M_gas)*T^1.75} / &
3794!      {p_atm/p_stand * (d_air^(1/3) + d_gas^(1/3))^2 }
3795! M_air = 28.965 : molar mass of air (g/mol)
3796! d_air = 19.70  : diffusion volume of air
3797! M_h2so4 = 98.08 : molar mass of h2so4 (g/mol)
3798! d_h2so4 = 51.96  : diffusion volume of h2so4
3799!
3800!> Called from main aerosol model
3801!
3802!> fxm: calculated for empty bins too
3803!> fxm: same diffusion coefficients and mean free paths used for sulphuric acid
3804!>      and organic vapours (average values? 'real' values for each?)
3805!> fxm: one should really couple with vapour production and loss terms as well
3806!>      should nucleation be coupled here as well????
3807!
3808! Coded by:
3809! Hannele Korhonen (FMI) 2005
3810! Harri Kokkola (FMI) 2006
3811! Juha Tonttila (FMI) 2014
3812! Rewritten to PALM by Mona Kurppa (UHel) 2017
3813!------------------------------------------------------------------------------!
3814 SUBROUTINE condensation( paero, pcsa, pcocnv, pcocsv, pchno3, pcnh3, pcw, pcs,&
3815                          ptemp, ppres, ptstep, prtcl )
3816       
3817    IMPLICIT NONE
3818   
3819!-- Input and output variables
3820    REAL(wp), INTENT(IN) ::  ppres !< ambient pressure (Pa)
3821    REAL(wp), INTENT(IN) ::  pcs   !< Water vapour saturation concentration
3822                                   !< (kg/m3)     
3823    REAL(wp), INTENT(IN) ::  ptemp !< ambient temperature (K)
3824    REAL(wp), INTENT(IN) ::  ptstep            !< timestep (s) 
3825    TYPE(component_index), INTENT(in) :: prtcl !< Keeps track which substances
3826                                               !< are used                                               
3827    REAL(wp), INTENT(INOUT) ::  pchno3 !< Gas concentrations (#/m3):
3828                                       !< nitric acid HNO3
3829    REAL(wp), INTENT(INOUT) ::  pcnh3  !< ammonia NH3
3830    REAL(wp), INTENT(INOUT) ::  pcocnv !< non-volatile organics
3831    REAL(wp), INTENT(INOUT) ::  pcocsv !< semi-volatile organics
3832    REAL(wp), INTENT(INOUT) ::  pcsa   !< sulphuric acid H2SO4
3833    REAL(wp), INTENT(INOUT) ::  pcw    !< Water vapor concentration (kg/m3)
3834    TYPE(t_section), INTENT(inout) ::  paero(fn2b) !< Aerosol properties                                     
3835!-- Local variables
3836    REAL(wp) ::  zbeta(fn2b) !< transitional correction factor for aerosols
3837    REAL(wp) ::  zcolrate(fn2b) !< collision rate of molecules to particles
3838                                !< (1/s)
3839    REAL(wp) ::  zcolrate_ocnv(fn2b) !< collision rate of organic molecules
3840                                     !< to particles (1/s)
3841    REAL(wp) ::  zcs_ocnv !< condensation sink of nonvolatile organics (1/s)       
3842    REAL(wp) ::  zcs_ocsv !< condensation sink of semivolatile organics (1/s)
3843    REAL(wp) ::  zcs_su !< condensation sink of sulfate (1/s)
3844    REAL(wp) ::  zcs_tot!< total condensation sink (1/s) (gases)
3845!-- vapour concentration after time step (#/m3)
3846    REAL(wp) ::  zcvap_new1 !< sulphuric acid
3847    REAL(wp) ::  zcvap_new2 !< nonvolatile organics
3848    REAL(wp) ::  zcvap_new3 !< semivolatile organics
3849    REAL(wp) ::  zdfpart(in1a+1) !< particle diffusion coefficient (m2/s)     
3850    REAL(wp) ::  zdfvap !< air diffusion coefficient (m2/s)
3851!-- change in vapour concentration (#/m3)
3852    REAL(wp) ::  zdvap1 !< sulphuric acid
3853    REAL(wp) ::  zdvap2 !< nonvolatile organics
3854    REAL(wp) ::  zdvap3 !< semivolatile organics
3855    REAL(wp) ::  zdvoloc(fn2b) !< change of organics volume in each bin [fxm]   
3856    REAL(wp) ::  zdvolsa(fn2b) !< change of sulphate volume in each bin [fxm]
3857    REAL(wp) ::  zj3n3(2)      !< Formation massrate of molecules in
3858                               !< nucleation, (molec/m3s). 1: H2SO4
3859                               !< and 2: organic vapor       
3860    REAL(wp) ::  zknud(fn2b) !< particle Knudsen number       
3861    REAL(wp) ::  zmfp    !< mean free path of condensing vapour (m)
3862    REAL(wp) ::  zrh     !< Relative humidity [0-1]         
3863    REAL(wp) ::  zvisc   !< viscosity of air (kg/(m s))     
3864    REAL(wp) ::  zn_vs_c !< ratio of nucleation of all mass transfer in the
3865                         !< smallest bin
3866    REAL(wp) ::  zxocnv  !< ratio of organic vapour in 3nm particles
3867    REAL(wp) ::  zxsa    !< Ratio in 3nm particles: sulphuric acid
3868   
3869    zj3n3  = 0.0_wp
3870    zrh    = pcw / pcs   
3871    zxocnv = 0.0_wp
3872    zxsa   = 0.0_wp
3873!
3874!-- Nucleation
3875    IF ( nsnucl > 0 )  THEN
3876       CALL nucleation( paero, ptemp, zrh, ppres, pcsa, pcocnv, pcnh3, ptstep, &
3877                        zj3n3, zxsa, zxocnv )
3878    ENDIF
3879!
3880!-- Condensation on pre-existing particles
3881    IF ( lscndgas )  THEN
3882!
3883!--    Initialise:
3884       zdvolsa = 0.0_wp 
3885       zdvoloc = 0.0_wp
3886       zcolrate = 0.0_wp
3887!             
3888!--    1) Properties of air and condensing gases:
3889!--    Viscosity of air (kg/(m s)) (Eq. 4.54 in Jabonson (2005))
3890       zvisc = ( 7.44523E-3_wp * ptemp ** 1.5_wp ) / ( 5093.0_wp *             &
3891                 ( ptemp + 110.4_wp ) )
3892!--    Diffusion coefficient of air (m2/s)
3893       zdfvap = 5.1111E-10_wp * ptemp ** 1.75_wp * ( p_0 + 1325.0_wp ) / ppres 
3894!--    Mean free path (m): same for H2SO4 and organic compounds
3895       zmfp = 3.0_wp * zdfvap * SQRT( pi * amh2so4 / ( 8.0_wp * argas * ptemp ) )
3896!                   
3897!--    2) Transition regime correction factor zbeta for particles:
3898!--       Fuchs and Sutugin (1971), In: Hidy et al. (ed.) Topics in current
3899!--       aerosol research, Pergamon. Size of condensing molecule considered 
3900!--       only for nucleation mode (3 - 20 nm)
3901!
3902!--    Particle Knudsen number: condensation of gases on aerosols
3903       zknud(in1a:in1a+1) = 2.0_wp * zmfp / ( paero(in1a:in1a+1)%dwet + d_sa )
3904       zknud(in1a+2:fn2b) = 2.0_wp * zmfp / paero(in1a+2:fn2b)%dwet
3905!   
3906!--    Transitional correction factor: aerosol + gas (the semi-empirical Fuchs-
3907!--    Sutugin interpolation function (Fuchs and Sutugin, 1971))
3908       zbeta = ( zknud + 1.0_wp ) / ( 0.377_wp * zknud + 1.0_wp + 4.0_wp /     &
3909               ( 3.0_wp * massacc ) * ( zknud + zknud ** 2.0_wp ) )
3910!                   
3911!--    3) Collision rate of molecules to particles
3912!--       Particle diffusion coefficient considered only for nucleation mode
3913!--       (3 - 20 nm)
3914!
3915!--    Particle diffusion coefficient (m2/s) (e.g. Eq. 15.29 in Jacobson (2005))
3916       zdfpart = abo * ptemp * zbeta(in1a:in1a+1) / ( 3.0_wp * pi * zvisc *    &
3917                 paero(in1a:in1a+1)%dwet )
3918!             
3919!--    Collision rate (mass-transfer coefficient): gases on aerosols (1/s)
3920!--    (Eq. 16.64 in Jacobson (2005))
3921       zcolrate(in1a:in1a+1) = MERGE( 2.0_wp * pi *                            &
3922                                      ( paero(in1a:in1a+1)%dwet + d_sa ) *     &
3923                                      ( zdfvap + zdfpart ) * zbeta(in1a:in1a+1)& 
3924                                        * paero(in1a:in1a+1)%numc, 0.0_wp,     &
3925                                      paero(in1a:in1a+1)%numc > nclim )
3926       zcolrate(in1a+2:fn2b) = MERGE( 2.0_wp * pi * paero(in1a+2:fn2b)%dwet *  &
3927                                      zdfvap * zbeta(in1a+2:fn2b) *            &
3928                                      paero(in1a+2:fn2b)%numc, 0.0_wp,         &
3929                                      paero(in1a+2:fn2b)%numc > nclim )
3930!                 
3931!-- 4) Condensation sink (1/s)
3932       zcs_tot = SUM( zcolrate )   ! total sink
3933!
3934!--    5) Changes in gas-phase concentrations and particle volume
3935!
3936!--    5.1) Organic vapours
3937!
3938!--    5.1.1) Non-volatile organic compound: condenses onto all bins
3939       IF ( pcocnv > 1.0E+10_wp  .AND.  zcs_tot > 1.0E-30_wp  .AND.            &
3940            is_used( prtcl,'OC' ) )                                            &
3941       THEN
3942!--       Ratio of nucleation vs. condensation rates in the smallest bin   
3943          zn_vs_c = 0.0_wp 
3944          IF ( zj3n3(2) > 1.0_wp )  THEN
3945             zn_vs_c = ( zj3n3(2) ) / ( zj3n3(2) + pcocnv * zcolrate(in1a) )
3946          ENDIF
3947!       
3948!--       Collision rate in the smallest bin, including nucleation and
3949!--       condensation(see Jacobson, Fundamentals of Atmospheric Modeling, 2nd
3950!--       Edition (2005), equation (16.73) )
3951          zcolrate_ocnv = zcolrate
3952          zcolrate_ocnv(in1a) = zcolrate_ocnv(in1a) + zj3n3(2) / pcocnv
3953!       
3954!--       Total sink for organic vapor
3955          zcs_ocnv = zcs_tot + zj3n3(2) / pcocnv
3956!       
3957!--       New gas phase concentration (#/m3)
3958          zcvap_new2 = pcocnv / ( 1.0_wp + ptstep * zcs_ocnv )
3959!       
3960!--       Change in gas concentration (#/m3)
3961          zdvap2 = pcocnv - zcvap_new2
3962!
3963!--       Updated vapour concentration (#/m3)               
3964          pcocnv = zcvap_new2
3965!       
3966!--       Volume change of particles (m3(OC)/m3(air))
3967          zdvoloc = zcolrate_ocnv(in1a:fn2b) / zcs_ocnv * amvoc * zdvap2
3968!       
3969!--       Change of volume due to condensation in 1a-2b
3970          paero(in1a:fn2b)%volc(2) = paero(in1a:fn2b)%volc(2) + zdvoloc 
3971!       
3972!--       Change of number concentration in the smallest bin caused by
3973!--       nucleation (Jacobson (2005), equation (16.75)). If zxocnv = 0, then 
3974!--       the chosen nucleation mechanism doesn't take into account the non-
3975!--       volatile organic vapors and thus the paero doesn't have to be updated.
3976          IF ( zxocnv > 0.0_wp )  THEN
3977             paero(in1a)%numc = paero(in1a)%numc + zn_vs_c * zdvoloc(in1a) /   &
3978                                amvoc / ( n3 * zxocnv )
3979          ENDIF
3980       ENDIF
3981!   
3982!--    5.1.2) Semivolatile organic compound: all bins except subrange 1
3983       zcs_ocsv = SUM( zcolrate(in2a:fn2b) ) !< sink for semi-volatile organics
3984       IF ( pcocsv > 1.0E+10_wp  .AND.  zcs_ocsv > 1.0E-30  .AND.              &
3985            is_used( prtcl,'OC') )                                             &
3986       THEN
3987!
3988!--       New gas phase concentration (#/m3)
3989          zcvap_new3 = pcocsv / ( 1.0_wp + ptstep * zcs_ocsv )
3990!       
3991!--       Change in gas concentration (#/m3)
3992          zdvap3 = pcocsv - zcvap_new3 
3993!       
3994!--       Updated gas co