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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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