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

Last change on this file since 3481 was 3481, checked in by raasch, 3 years ago

temporary variable cc introduced to circumvent a possible Intel18 compiler bug related to contiguous/non-contguous pointer/target attributes

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