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

Last change on this file since 3767 was 3767, checked in by raasch, 2 years ago

unused variables removed from rrd-subroutines parameter list

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