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

Last change on this file since 3833 was 3833, checked in by forkel, 2 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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