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

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

Bugfix for previous commit

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