Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 3956)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4012)
@@ -26,4 +26,16 @@
 ! -----------------
 ! $Id$
+! Merge salsa branch to trunk. List of changes:
+! - Error corrected in distr_update that resulted in the aerosol number size
+!   distribution not converging if the concentration was nclim.
+! - Added a separate output for aerosol liquid water (s_H2O)
+! - aerosol processes for a size bin are now calculated only if the aerosol
+!   number of concentration of that bin is > 2*nclim
+! - An initialisation error in the subroutine "deposition" corrected and the
+!   subroutine reformatted.
+! - stuff from salsa_util_mod.f90 moved into salsa_mod.f90
+! - calls for closing the netcdf input files added
+! 
+! 3956 2019-05-07 12:32:52Z monakurppa
 ! - Conceptual bug in depo_surf correct for urban and land surface model
 ! - Subroutine salsa_tendency_ij optimized.
@@ -147,30 +159,30 @@
 !> @todo emission mode "parameterized", i.e. based on street type
 !> @todo Allow insoluble emissions
-!> @todo two-way nesting is not working properly
+!> @todo Apply flux limiter in prognostic equations
 !------------------------------------------------------------------------------!
  MODULE salsa_mod
 
-    USE basic_constants_and_equations_mod,                                     &
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  c_p, g, p_0, pi, r_d
 
-    USE chem_gasphase_mod,                                                     &
+    USE chem_gasphase_mod,                                                                         &
         ONLY:  nspec, nvar, spc_names
 
-    USE chem_modules,                                                          &
+    USE chem_modules,                                                                              &
         ONLY:  call_chem_at_all_substeps, chem_gasphase_on, chem_species
 
     USE control_parameters
 
-    USE indices,                                                               &
-        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
-               nzb_s_inner, nz, nzt, wall_flags_0
+    USE indices,                                                                                   &
+        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nz, nzt, wall_flags_0
 
     USE kinds
 
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  chem_emis_att_type, chem_emis_val_type
+
     USE pegrid
 
-    USE salsa_util_mod
-
-    USE statistics,                                                            &
+    USE statistics,                                                                                &
         ONLY:  sums_salsa_ws_l
 
@@ -452,5 +464,5 @@
                                                                (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
     REAL(wp), DIMENSION(nmod) ::  surface_aerosol_flux = &
-                             (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
+                                 (/1.0E+8_wp, 1.0E+9_wp, 1.0E+5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  bin_low_limits     !< to deliver information about
@@ -464,4 +476,11 @@
 !
 !-- SALSA derived datatypes:
+!
+!-- Component index
+    TYPE component_index
+       CHARACTER(len=3), ALLOCATABLE ::  comp(:)  !< Component name
+       INTEGER(iwp) ::  ncomp  !< Number of components
+       INTEGER(iwp), ALLOCATABLE ::  ind(:)  !< Component index
+    END TYPE component_index
 !
 !-- For matching LSM and USM surface types and the deposition module surface types
@@ -607,4 +626,7 @@
     TYPE(salsa_emission_value_type)     ::  aero_emission      !< emission values
     TYPE(salsa_emission_mode_type)      ::  def_modes          !< default emission modes
+
+    TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
+    TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
 
     TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
@@ -979,4 +1001,7 @@
  SUBROUTINE salsa_header( io )
 
+    USE indices,                                                                                   &
+        ONLY:  nx, ny, nz
+
     IMPLICIT NONE
  
@@ -987,7 +1012,7 @@
     WRITE( io, 2 ) skip_time_do_salsa
     WRITE( io, 3 ) dt_salsa
-    WRITE( io, 4 )  SHAPE( aerosol_number(1)%conc ), nbins_aerosol
+    WRITE( io, 4 )  nz, ny, nx, nbins_aerosol
     IF ( advect_particle_water )  THEN
-       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol,             &
+       WRITE( io, 5 )  SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol,              &
                         advect_particle_water
     ELSE
@@ -1154,4 +1179,5 @@
                  aerosol_number(i)%init(nzb:nzt+1),                                                &
                  aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+       aerosol_number(i)%init = nclim
        IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
           ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
@@ -1179,4 +1205,5 @@
                  aerosol_mass(i)%init(nzb:nzt+1),                                                  &
                  aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
+       aerosol_mass(i)%init = mclim
        IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
           ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
@@ -1283,4 +1310,5 @@
                     salsa_gas(i)%init(nzb:nzt+1),                              &
                     salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+          salsa_gas(i)%init = nclim
           IF ( include_emission )  ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
        ENDDO
@@ -1402,28 +1430,12 @@
     IF ( nldepo )  sedim_vd = 0.0_wp
 
-    DO  ib = 1, nbins_aerosol
-       IF ( .NOT. read_restart_data_salsa )  aerosol_number(ib)%conc = nclim
-       aerosol_number(ib)%conc_p    = 0.0_wp
-       aerosol_number(ib)%tconc_m   = 0.0_wp
-       aerosol_number(ib)%flux_s    = 0.0_wp
-       aerosol_number(ib)%diss_s    = 0.0_wp
-       aerosol_number(ib)%flux_l    = 0.0_wp
-       aerosol_number(ib)%diss_l    = 0.0_wp
-       aerosol_number(ib)%init      = nclim
-       aerosol_number(ib)%sums_ws_l = 0.0_wp
-    ENDDO
-    DO  ic = 1, ncomponents_mass*nbins_aerosol
-       IF ( .NOT. read_restart_data_salsa )  aerosol_mass(ic)%conc = mclim
-       aerosol_mass(ic)%conc_p    = 0.0_wp
-       aerosol_mass(ic)%tconc_m   = 0.0_wp
-       aerosol_mass(ic)%flux_s    = 0.0_wp
-       aerosol_mass(ic)%diss_s    = 0.0_wp
-       aerosol_mass(ic)%flux_l    = 0.0_wp
-       aerosol_mass(ic)%diss_l    = 0.0_wp
-       aerosol_mass(ic)%init      = mclim
-       aerosol_mass(ic)%sums_ws_l = 0.0_wp
-    ENDDO
-
     IF ( .NOT. salsa_gases_from_chem )  THEN
+       IF ( .NOT. read_restart_data_salsa )  THEN
+          salsa_gas(1)%conc = h2so4_init
+          salsa_gas(2)%conc = hno3_init
+          salsa_gas(3)%conc = nh3_init
+          salsa_gas(4)%conc = ocnv_init
+          salsa_gas(5)%conc = ocsv_init
+       ENDIF
        DO  ig = 1, ngases_salsa
           salsa_gas(ig)%conc_p    = 0.0_wp
@@ -1434,14 +1446,8 @@
           salsa_gas(ig)%diss_l    = 0.0_wp
           salsa_gas(ig)%sums_ws_l = 0.0_wp
+          salsa_gas(ig)%conc_p    = salsa_gas(ig)%conc
        ENDDO
-       IF ( .NOT. read_restart_data_salsa )  THEN
-          salsa_gas(1)%conc = h2so4_init
-          salsa_gas(2)%conc = hno3_init
-          salsa_gas(3)%conc = nh3_init
-          salsa_gas(4)%conc = ocnv_init
-          salsa_gas(5)%conc = ocsv_init 
-       ENDIF
-!
-!--    Set initial value for gas compound tracers and initial values
+!
+!--    Set initial value for gas compound tracer
        salsa_gas(1)%init = h2so4_init
        salsa_gas(2)%init = hno3_init
@@ -1466,5 +1472,5 @@
 !-- Initialise location-dependent aerosol size distributions and chemical compositions:
     CALL aerosol_init
-!
+
 !-- Initalisation run of SALSA + calculate the vertical top index of the topography
     DO  i = nxl, nxr
@@ -1477,4 +1483,24 @@
        ENDDO
     ENDDO
+
+    DO  ib = 1, nbins_aerosol
+       aerosol_number(ib)%conc_p    = aerosol_number(ib)%conc
+       aerosol_number(ib)%tconc_m   = 0.0_wp
+       aerosol_number(ib)%flux_s    = 0.0_wp
+       aerosol_number(ib)%diss_s    = 0.0_wp
+       aerosol_number(ib)%flux_l    = 0.0_wp
+       aerosol_number(ib)%diss_l    = 0.0_wp
+       aerosol_number(ib)%sums_ws_l = 0.0_wp
+    ENDDO
+    DO  ic = 1, ncomponents_mass*nbins_aerosol
+       aerosol_mass(ic)%conc_p    = aerosol_mass(ic)%conc
+       aerosol_mass(ic)%tconc_m   = 0.0_wp
+       aerosol_mass(ic)%flux_s    = 0.0_wp
+       aerosol_mass(ic)%diss_s    = 0.0_wp
+       aerosol_mass(ic)%flux_l    = 0.0_wp
+       aerosol_mass(ic)%diss_l    = 0.0_wp
+       aerosol_mass(ic)%sums_ws_l = 0.0_wp
+    ENDDO
+!
 !
 !-- Initialise the deposition scheme and surface types
@@ -1587,5 +1613,6 @@
 
     USE netcdf_data_input_mod,                                                                     &
-        ONLY:  get_attribute, get_variable, netcdf_data_input_get_dimension_length, open_read_file
+        ONLY:  close_input_file, get_attribute, get_variable,                                      &
+               netcdf_data_input_get_dimension_length, open_read_file
 
     IMPLICIT NONE
@@ -1615,5 +1642,4 @@
 
     REAL(wp), DIMENSION(nbins_aerosol) ::  core   !< size of the bin mid aerosol particle
-    REAL(wp), DIMENSION(nbins_aerosol) ::  nsect  !< size distribution (#/m3)
 
     REAL(wp), DIMENSION(0:nz+1) ::  pnf2a   !< number fraction in 2a
@@ -1638,5 +1664,4 @@
 !
 !-- Set concentrations to zero
-    nsect(:)     = 0.0_wp
     pndist(:,:)  = 0.0_wp
     pnf2a(:)     = nf2a
@@ -1667,5 +1692,5 @@
 !
 !--       Allocate memory
-          ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                            &
+          ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc),                              &
                     pr_mass_fracs_b(nzb:nzt+1,pr_ncc) )
           pr_mass_fracs_a = 0.0_wp
@@ -1806,5 +1831,7 @@
                            ' for SALSA missing!'
           CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 )
-
+!
+!--       Close input file
+          CALL close_input_file( id_dyn )
        ENDIF   ! netcdf_extend
 
@@ -1892,7 +1919,9 @@
              salsa_gas(ig)%init(nzb)   =  salsa_gas(ig)%init(nzb+1)
              salsa_gas(ig)%init(nzt+1) =  salsa_gas(ig)%init(nzt)
-             DO  k = nzb, nzt+1
-                salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
-             ENDDO
+             IF ( .NOT. read_restart_data_salsa )  THEN
+                DO  k = nzb, nzt+1
+                   salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k)
+                ENDDO
+             ENDIF
           ENDDO
 
@@ -1901,4 +1930,7 @@
                            ' for SALSA missing!'
           CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 )
+!
+!--       Close input file
+          CALL close_input_file( id_dyn )
        ENDIF   ! netcdf_extend
 #else
@@ -1940,5 +1972,7 @@
 !--          Region 1:
              DO  ib = start_subrange_1a, end_subrange_1a
-                aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
+                IF ( .NOT. read_restart_data_salsa )  THEN
+                   aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag
+                ENDIF
                 IF ( prunmode == 1 )  THEN
                    aerosol_number(ib)%init = pndist(:,ib)
@@ -1949,5 +1983,7 @@
              IF ( nreg > 1 )  THEN
                 DO  ib = start_subrange_2a, end_subrange_2a
-                   aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
+                   IF ( .NOT. read_restart_data_salsa )  THEN
+                      aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag
+                   ENDIF
                    IF ( prunmode == 1 )  THEN
                       aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib)
@@ -1957,6 +1993,8 @@
                    DO  ib = start_subrange_2b, end_subrange_2b
                       IF ( pnf2a(k) < 1.0_wp )  THEN
-                         aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *       &
-                                                          pndist(k,ib) * flag
+                         IF ( .NOT. read_restart_data_salsa )  THEN
+                            aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pnf2a(k) ) *    &
+                                                             pndist(k,ib) * flag
+                         ENDIF
                          IF ( prunmode == 1 )  THEN
                             aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp - nf2a ) * pndist(:,ib)
@@ -1976,6 +2014,8 @@
                 ib = start_subrange_1a
                 DO  ic = ss, ee
-                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib)&
-                                                  * core(ib) * arhoh2so4 * flag
+                   IF ( .NOT. read_restart_data_salsa )  THEN
+                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) *          &
+                                                     pndist(k,ib) * core(ib) * arhoh2so4 * flag
+                   ENDIF
                    IF ( prunmode == 1 )  THEN
                       aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp - pmfoc1a(k) ) * pndist(k,ib) &
@@ -1991,7 +2031,9 @@
                 ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a !< end
                 ib = start_subrange_1a
-                DO  ic = ss, ee 
-                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *       &
-                                                  core(ib) * arhooc * flag
+                DO  ic = ss, ee
+                   IF ( .NOT. read_restart_data_salsa )  THEN
+                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *    &
+                                                     core(ib) * arhooc * flag
+                   ENDIF
                    IF ( prunmode == 1 )  THEN
                       aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) *        &
@@ -2143,6 +2185,8 @@
              ib = start_subrange_2a
              DO ic = ss, ee
-                aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) * &
-                                              pcore(ib) * prho * flag
+                IF ( .NOT. read_restart_data_salsa )  THEN
+                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib)&
+                                                  * pcore(ib) * prho * flag
+                ENDIF
                 IF ( prunmode == 1 )  THEN
                    aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) *  &
@@ -2158,6 +2202,8 @@
                 ib = start_subrange_2a
                 DO ic = ss, ee
-                   aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) *&
-                                                  pndist(k,ib) * pcore(ib) * prho * flag
+                   IF ( .NOT. read_restart_data_salsa )  THEN
+                      aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k))&
+                                                     * pndist(k,ib) * pcore(ib) * prho * flag
+                   ENDIF
                    IF ( prunmode == 1 )  THEN
                       aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp - pnf2a(k) ) * &
@@ -2810,7 +2856,6 @@
        ENDDO
 !
-!--    On EACH call of salsa_driver, calculate the ambient sizes of 
-!--    particles by equilibrating soluble fraction of particles with water
-!--    using the ZSR method.
+!--    Calculate the ambient sizes of particles by equilibrating soluble fraction of particles with
+!--    water using the ZSR method.
        in_rh = in_cw(k) / in_cs(k)
        IF ( prunmode==1  .OR.  .NOT. advect_particle_water )  THEN
@@ -2865,5 +2910,5 @@
 !--    Calculate changes in concentrations
        DO ib = 1, nbins_aerosol
-          aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc -      &
+          aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( lo_aero(ib)%numc -   &
                                            aero_old(ib)%numc ) * flag
        ENDDO
@@ -2875,5 +2920,5 @@
           ic = 1
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhoh2so4 * flag
              ic = ic+1
@@ -2887,5 +2932,5 @@
           ic = 1
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhooc * flag
              ic = ic+1
@@ -2899,5 +2944,5 @@
           ic = 1 + end_subrange_1a
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhobc * flag
              ic = ic+1
@@ -2911,5 +2956,5 @@
           ic = 1 + end_subrange_1a
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhodu * flag
              ic = ic+1
@@ -2923,5 +2968,5 @@
           ic = 1 + end_subrange_1a
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhoss * flag
              ic = ic+1
@@ -2935,5 +2980,5 @@
           ic = 1
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhohno3 * flag
              ic = ic+1
@@ -2947,5 +2992,5 @@
           ic = 1
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhonh3 * flag
              ic = ic+1
@@ -2960,14 +3005,23 @@
           ic = 1
           DO ss = str, endi
-             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -   &
+             aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( lo_aero(ic)%volc(vc) -&
                                             aero_old(ic)%volc(vc) ) * arhoh2o * flag
-             IF ( prunmode == 1 )  THEN
-                aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k),                          &
-                                                aerosol_mass(ss)%conc(k,j,i) )
-                IF ( k == nzb+1 )  THEN
-                   aerosol_mass(ss)%init(k-1) = 0.0_wp
-                ELSEIF ( k == nzt  )  THEN
-                   aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
-                ENDIF
+             ic = ic+1
+          ENDDO
+       ENDIF
+       IF ( prunmode == 1 )  THEN
+          nc_h2o = get_index( prtcl,'H2O' )
+          vc = 8
+          str = ( nc_h2o-1 ) * nbins_aerosol + 1
+          endi = nc_h2o * nbins_aerosol
+          ic = 1
+          DO ss = str, endi
+             aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
+                                             aero_old(ic)%volc(vc) ) * arhoh2o )
+             IF ( k == nzb+1 )  THEN
+                aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
+             ELSEIF ( k == nzt  )  THEN
+                aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
+                aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
              ENDIF
              ic = ic+1
@@ -3446,5 +3500,4 @@
     REAL(wp) ::  avis              !< molecular viscocity of air (kg/(m*s))
     REAL(wp) ::  beta_im           !< parameter for turbulent impaction
-    REAL(wp) ::  beta              !< Cunningham slip-flow correction factor
     REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
     REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
@@ -3453,5 +3506,4 @@
     REAL(wp) ::  depo              !< deposition velocity (m/s)
     REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
-    REAL(wp) ::  Kn                !< Knudsen number
     REAL(wp) ::  lambda            !< molecular mean free path (m)
     REAL(wp) ::  mdiff             !< particle diffusivity coefficient
@@ -3461,6 +3513,5 @@
     REAL(wp) ::  pdn               !< particle density (kg/m3)
     REAL(wp) ::  ustar             !< friction velocity (m/s)
-    REAL(wp) ::  va                !< thermal speed of an air molecule (m/s) 
-    REAL(wp) ::  zdwet             !< wet diameter (m)
+    REAL(wp) ::  va                !< thermal speed of an air molecule (m/s)
 
     REAL(wp), INTENT(in) ::  adn    !< air density (kg/m3)
@@ -3471,4 +3522,8 @@
     REAL(wp), INTENT(inout) ::  kvis   !< kinematic viscosity of air (m2/s)
 
+    REAL(wp), DIMENSION(nbins_aerosol) ::  beta   !< Cunningham slip-flow correction factor
+    REAL(wp), DIMENSION(nbins_aerosol) ::  Kn     !< Knudsen number
+    REAL(wp), DIMENSION(nbins_aerosol) ::  zdwet  !< wet diameter (m)
+
     REAL(wp), DIMENSION(:), INTENT(inout) ::  schmidt_num  !< particle Schmidt number
     REAL(wp), DIMENSION(:), INTENT(inout) ::  vc  !< critical fall speed i.e. settling velocity of
@@ -3494,39 +3549,43 @@
     lambda = 2.0_wp * avis / ( adn * va )
 !
-!-- Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
-    alpha   = alpha_z01(depo_pcm_type_num)
-    gamma   = gamma_z01(depo_pcm_type_num)
-    par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
-!
-!-- Deposition efficiencies from Table 1. Constants from Table 2.
-    par_l            = l_p10(depo_pcm_type_num) * 0.01_wp
-    c_brownian_diff  = c_b_p10(depo_pcm_type_num)
-    c_interception   = c_in_p10(depo_pcm_type_num)
-    c_impaction      = c_im_p10(depo_pcm_type_num)
-    beta_im          = beta_im_p10(depo_pcm_type_num)
-    c_turb_impaction = c_it_p10(depo_pcm_type_num)
-
-    DO  ib = 1, nbins_aerosol
-
-       IF ( paero(ib)%numc < nclim )  CYCLE
-       zdwet = paero(ib)%dwet
-!
-!--    Knudsen number (Eq. 15.23)
-       Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
-!
-!--    Cunningham slip-flow correction (Eq. 15.30)
-       beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
-
-!--    Particle diffusivity coefficient (Eq. 15.29)
-       mdiff = ( abo * tk * beta ) / ( 3.0_wp * pi * avis * zdwet )
-!
-!--    Particle Schmidt number (Eq. 15.36)
-       schmidt_num(ib) = kvis / mdiff
-!
-!--    Critical fall speed i.e. settling velocity  (Eq. 20.4)
-       vc(ib) = MIN( 1.0_wp, terminal_vel( 0.5_wp * zdwet, pdn, adn, avis, beta) )
-!
-!--    Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
-       IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
+!-- Particle wet diameter (m)
+    zdwet = paero(:)%dwet
+!
+!-- Knudsen number (Eq. 15.23)
+    Kn = MAX( 1.0E-2_wp, lambda / ( zdwet * 0.5_wp ) ) ! To avoid underflow
+!
+!-- Cunningham slip-flow correction (Eq. 15.30)
+    beta = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
+!
+!-- Critical fall speed i.e. settling velocity  (Eq. 20.4)
+    vc = MIN( 1.0_wp, zdwet**2 * ( pdn - adn ) * g * beta / ( 18.0_wp * avis ) )
+!
+!-- Deposition on vegetation
+    IF ( lsdepo_pcm  .AND.  plant_canopy  .AND.  lad > 0.0_wp )  THEN
+!
+!--    Parameters for the land use category 'deciduous broadleaf trees'(Table 3)
+       alpha   = alpha_z01(depo_pcm_type_num)
+       gamma   = gamma_z01(depo_pcm_type_num)
+       par_a   = A_z01(depo_pcm_type_num, season) * 1.0E-3_wp
+!
+!--    Deposition efficiencies from Table 1. Constants from Table 2.
+       par_l            = l_p10(depo_pcm_type_num) * 0.01_wp
+       c_brownian_diff  = c_b_p10(depo_pcm_type_num)
+       c_interception   = c_in_p10(depo_pcm_type_num)
+       c_impaction      = c_im_p10(depo_pcm_type_num)
+       beta_im          = beta_im_p10(depo_pcm_type_num)
+       c_turb_impaction = c_it_p10(depo_pcm_type_num)
+
+       DO  ib = 1, nbins_aerosol
+
+          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
+
+!--       Particle diffusivity coefficient (Eq. 15.29)
+          mdiff = ( abo * tk * beta(ib) ) / ( 3.0_wp * pi * avis * zdwet(ib) )
+!
+!--       Particle Schmidt number (Eq. 15.36)
+          schmidt_num(ib) = kvis / mdiff
+!
+!--       Friction velocity for deposition on vegetation. Calculated following Prandtl (1925):
           ustar = SQRT( cdc ) * mag_u
           SELECT CASE ( depo_pcm_par_num )
@@ -3546,6 +3605,7 @@
              paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa
           ENDDO
-       ENDIF
-    ENDDO
+       ENDDO
+
+    ENDIF
 
  END SUBROUTINE deposition
@@ -3779,6 +3839,6 @@
 
              DO  ib = 1, nbins_aerosol
-                IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
-                   CYCLE
+                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
+                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
 
                 SELECT CASE ( depo_surf_par_num )
@@ -3810,6 +3870,6 @@
 
              DO  ib = 1, nbins_aerosol
-                IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
-                   CYCLE
+                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
+                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
 
                 SELECT CASE ( depo_surf_par_num )
@@ -3841,6 +3901,6 @@
 
              DO  ib = 1, nbins_aerosol
-                IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
-                   CYCLE
+                IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                     &
+                     schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
 
                 SELECT CASE ( depo_surf_par_num )
@@ -3860,4 +3920,5 @@
 
           DO  ib = 1, nbins_aerosol
+             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim ) )  CYCLE
 !
 !--          Calculate changes in surface fluxes due to dry deposition
@@ -3891,6 +3952,6 @@
 
           DO  ib = 1, nbins_aerosol
-             IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp )    &
-                CYCLE
+             IF ( aerosol_number(ib)%conc(k,j,i) < ( 2.0_wp * nclim )  .OR.                        &
+                  schmidt_num(k+1,ib) < 1.0_wp )  CYCLE
 
              SELECT CASE ( depo_surf_par_num )
@@ -3932,25 +3993,4 @@
 ! Description:
 ! ------------
-! Function for calculating terminal velocities for different particles sizes.
-!------------------------------------------------------------------------------!
- REAL(wp) FUNCTION terminal_vel( radius, rhop, rhoa, visc, beta )
-
-    IMPLICIT NONE
-
-    REAL(wp), INTENT(in) ::  beta    !< Cunningham correction factor
-    REAL(wp), INTENT(in) ::  radius  !< particle radius (m)
-    REAL(wp), INTENT(in) ::  rhop    !< particle density (kg/m3)
-    REAL(wp), INTENT(in) ::  rhoa    !< air density (kg/m3)
-    REAL(wp), INTENT(in) ::  visc    !< molecular viscosity of air (kg/(m*s))
-
-!
-!-- Stokes law with Cunningham slip correction factor
-    terminal_vel = 4.0_wp * radius**2 * ( rhop - rhoa ) * g * beta / ( 18.0_wp * visc ) ! (m/s)
-
- END FUNCTION terminal_vel
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Calculates particle loss and change in size distribution due to (Brownian)
 !> coagulation. Only for particles with dwet < 30 micrometres. 
@@ -4034,7 +4074,7 @@
 !-- Aero-aero coagulation
     DO  mm = 1, end_subrange_2b   ! smaller colliding particle
-       IF ( paero(mm)%numc < nclim )  CYCLE
+       IF ( paero(mm)%numc < ( 2.0_wp * nclim ) )  CYCLE
        DO  nn = mm, end_subrange_2b   ! larger colliding particle
-          IF ( paero(nn)%numc < nclim )  CYCLE
+          IF ( paero(nn)%numc < ( 2.0_wp * nclim ) )  CYCLE
 
           zdpart_mm = MIN( paero(mm)%dwet, 30.0E-6_wp )     ! Limit to 30 um
@@ -4053,5 +4093,5 @@
 !-- Aerosols in subrange 1a:
     DO  ib = start_subrange_1a, end_subrange_1a
-       IF ( paero(ib)%numc < nclim )  CYCLE
+       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
        zminusterm   = 0.0_wp
        zplusterm(:) = 0.0_wp
@@ -4080,5 +4120,5 @@
 !-- Aerosols in subrange 2a:
     DO  ib = start_subrange_2a, end_subrange_2a
-       IF ( paero(ib)%numc < nclim )  CYCLE
+       IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
        zminusterm   = 0.0_wp
        zplusterm(:) = 0.0_wp
@@ -4128,5 +4168,5 @@
     IF ( .NOT. no_insoluble )  THEN
        DO  ib = start_subrange_2b, end_subrange_2b
-          IF ( paero(ib)%numc < nclim )  CYCLE
+          IF ( paero(ib)%numc < ( 2.0_wp * nclim ) )  CYCLE
           zminusterm   = 0.0_wp
           zplusterm(:) = 0.0_wp
@@ -4337,14 +4377,14 @@
     REAL(wp), INTENT(inout) ::  pcw      !< Water vapor concentration (kg/m3)
 
-    REAL(wp), DIMENSION(nbins_aerosol)               ::  zbeta          !< transitional correction factor
-    REAL(wp), DIMENSION(nbins_aerosol)               ::  zcolrate       !< collision rate (1/s)
-    REAL(wp), DIMENSION(nbins_aerosol)               ::  zcolrate_ocnv  !< collision rate of non-vol. OC (1/s)
-    REAL(wp), DIMENSION(start_subrange_1a+1) ::  zdfpart        !< particle diffusion coefficient (m2/s)
-    REAL(wp), DIMENSION(nbins_aerosol)               ::  zdvoloc        !< change of organics volume
-    REAL(wp), DIMENSION(nbins_aerosol)               ::  zdvolsa        !< change of sulphate volume
+    REAL(wp), DIMENSION(nbins_aerosol)       ::  zbeta          !< transitional correction factor
+    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate       !< collision rate (1/s)
+    REAL(wp), DIMENSION(nbins_aerosol)       ::  zcolrate_ocnv  !< collision rate of OCNV (1/s)
+    REAL(wp), DIMENSION(start_subrange_1a+1) ::  zdfpart        !< particle diffusion coef. (m2/s)
+    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvoloc        !< change of organics volume
+    REAL(wp), DIMENSION(nbins_aerosol)       ::  zdvolsa        !< change of sulphate volume
     REAL(wp), DIMENSION(2)                   ::  zj3n3          !< Formation massrate of molecules
                                                                 !< in nucleation, (molec/m3s),
                                                                 !< 1: H2SO4 and 2: organic vapor
-    REAL(wp), DIMENSION(nbins_aerosol)   ::  zknud          !< particle Knudsen number
+    REAL(wp), DIMENSION(nbins_aerosol)       ::  zknud          !< particle Knudsen number
 
     TYPE(component_index), INTENT(in) :: prtcl  !< Keeps track which substances are used
@@ -7095,9 +7135,11 @@
     REAL(wp) ::  znfrac  !< number fraction to be moved to the larger bin
     REAL(wp) ::  zvfrac  !< volume fraction to be moved to the larger bin
-    REAL(wp) ::  zVexc   !< Volume in the grown bin which exceeds the bin upper limit
-    REAL(wp) ::  zVihi   !< particle volume at the high end of the bin
-    REAL(wp) ::  zVilo   !< particle volume at the low end of the bin
+    REAL(wp) ::  zvexc   !< Volume in the grown bin which exceeds the bin upper limit
+    REAL(wp) ::  zvihi   !< particle volume at the high end of the bin
+    REAL(wp) ::  zvilo   !< particle volume at the low end of the bin
     REAL(wp) ::  zvpart  !< particle volume (m3)
-    REAL(wp) ::  zVrat   !< volume ratio of a size bin
+    REAL(wp) ::  zvrat   !< volume ratio of a size bin
+
+    real(wp), dimension(nbins_aerosol) ::  dummy
 
     TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero !< aerosol properties
@@ -7106,4 +7148,6 @@
     zvfrac      = 0.0_wp
     within_bins = .FALSE.
+
+    dummy = paero(:)%numc
 !
 !-- Check if the volume of the bin is within bin limits after update
@@ -7150,20 +7194,20 @@
 !
 !--          Volume ratio of the size bin
-             zVrat = paero(ib)%vhilim / paero(ib)%vlolim
+             zvrat = paero(ib)%vhilim / paero(ib)%vlolim
 !
 !--          Particle volume at the low end of the bin
-             zVilo = 2.0_wp * zvpart / ( 1.0_wp + zVrat )
+             zvilo = 2.0_wp * zvpart / ( 1.0_wp + zvrat )
 !
 !--          Particle volume at the high end of the bin
-             zVihi = zVrat * zVilo
+             zvihi = zvrat * zvilo
 !
 !--          Volume in the grown bin which exceeds the bin upper limit
-             zVexc = 0.5_wp * ( zVihi + paero(ib)%vhilim )
+             zvexc = 0.5_wp * ( zvihi + paero(ib)%vhilim )
 !
 !--          Number fraction to be moved to the larger bin
-             znfrac = MIN( 1.0_wp, ( zVihi - paero(ib)%vhilim) / ( zVihi - zVilo ) )
+             znfrac = MIN( 1.0_wp, ( zvihi - paero(ib)%vhilim) / ( zvihi - zvilo ) )
 !
 !--          Volume fraction to be moved to the larger bin
-             zvfrac = MIN( 0.99_wp, znfrac * zVexc / zvpart )
+             zvfrac = MIN( 0.99_wp, znfrac * zvexc / zvpart )
              IF ( zvfrac < 0.0_wp )  THEN
                 message_string = 'Error: zvfrac < 0'
@@ -7175,8 +7219,8 @@
 !
 !--          Volume (cm3/cm3)
-             paero(mm)%volc(:) = paero(mm)%volc(:) + znfrac * paero(ib)%numc * zVexc *             &
-                                 paero(ib)%volc(:) / SUM( paero(ib)%volc(:) )
-             paero(ib)%volc(:) = paero(ib)%volc(:) - znfrac * paero(ib)%numc * zVexc *             &
-                                 paero(ib)%volc(:) / SUM( paero(ib)%volc(:) )
+             paero(mm)%volc(:) = paero(mm)%volc(:) + znfrac * paero(ib)%numc * zvexc *             &
+                                 paero(ib)%volc(:) / SUM( paero(ib)%volc(1:7) )
+             paero(ib)%volc(:) = paero(ib)%volc(:) - znfrac * paero(ib)%numc * zvexc *             &
+                                 paero(ib)%volc(:) / SUM( paero(ib)%volc(1:7) )
 
 !--          Number concentration (#/m3)
@@ -7266,11 +7310,12 @@
 !--    Calculate total mass concentration per bin
        mcsum = 0.0_wp
-       DO  ic = 1, ncc
+       DO  ic = 1, ncomponents_mass
           icc = ( ic - 1 ) * nbins_aerosol + ib
           mcsum = mcsum + aerosol_mass(icc)%conc(:,j,i) * flag
+          aerosol_mass(icc)%conc(:,j,i) = MAX( mclim, aerosol_mass(icc)%conc(:,j,i) ) * flag
        ENDDO
 !
 !--    Check that number and mass concentration match qualitatively
-       IF ( ANY ( aerosol_number(ib)%conc(:,j,i) >= nclim  .AND. mcsum <= 0.0_wp ) )  THEN
+       IF ( ANY( aerosol_number(ib)%conc(:,j,i) > nclim  .AND. mcsum <= 0.0_wp ) )  THEN
           DO  k = nzb+1, nzt
              IF ( aerosol_number(ib)%conc(k,j,i) >= nclim  .AND. mcsum(k) <= 0.0_wp )  THEN
@@ -7349,5 +7394,5 @@
                                                  flag_zddry > 0.0_wp )
        ENDDO
-       aerosol_number(ib)%conc(:,j,i) = MERGE( nclim, aerosol_number(ib)%conc(:,j,i),              &
+       aerosol_number(ib)%conc(:,j,i) = MERGE( nclim * flag, aerosol_number(ib)%conc(:,j,i),       &
                                                flag_zddry > 0.0_wp )
        ra_dry(:,j,i,ib) = MAX( 1.0E-10_wp, 0.5_wp * zddry )
@@ -8387,6 +8432,7 @@
 
     USE netcdf_data_input_mod,                                                                     &
-        ONLY:  check_existence, get_attribute, get_variable, inquire_num_variables,                &
-               inquire_variable_names, netcdf_data_input_get_dimension_length, open_read_file
+        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
+               inquire_num_variables, inquire_variable_names,                                      &
+               netcdf_data_input_get_dimension_length, open_read_file
 
     USE surface_mod,                                                                               &
@@ -8482,17 +8528,17 @@
 
 #if defined( __netcdf )
+!
+!--       Check existence of PIDS_SALSA file
+          INQUIRE( FILE = TRIM( input_file_salsa ) // TRIM( coupling_char ), EXIST = netcdf_extend )
+          IF ( .NOT. netcdf_extend )  THEN
+             message_string = 'Input file '// TRIM( input_file_salsa ) //  TRIM( coupling_char )&
+                              // ' missing!'
+             CALL message( 'salsa_emission_setup', 'PA0629', 1, 2, 0, 6, 0 )
+          ENDIF
+!
+!--       Open file in read-only mode
+          CALL open_read_file( TRIM( input_file_salsa ) // TRIM( coupling_char ), id_salsa )
+
           IF ( init )  THEN
-!
-!--          Check existence of PIDS_SALSA file
-             INQUIRE( FILE = TRIM( input_file_salsa ) // TRIM( coupling_char ),                    &
-                      EXIST = netcdf_extend )
-             IF ( .NOT. netcdf_extend )  THEN
-                message_string = 'Input file '// TRIM( input_file_salsa ) //  TRIM( coupling_char )&
-                                 // ' missing!'
-                CALL message( 'salsa_emission_setup', 'PA0629', 1, 2, 0, 6, 0 )
-             ENDIF
-!
-!--          Open file in read-only mode
-             CALL open_read_file( TRIM( input_file_salsa ) // TRIM( coupling_char ), id_salsa )
 !
 !--          Read the index and name of chemical components
@@ -8822,5 +8868,5 @@
 
 
-             DEALLOCATE( source_array )
+             DEALLOCATE( nsect_emission, source_array )
 !
 !--       Pre-processed:
@@ -8857,5 +8903,7 @@
 
           ENDIF
-
+!
+!--       Close input file
+          CALL close_input_file( id_salsa )
 #else
           message_string = 'salsa_emission_mode = "read_from_file", but preprocessor directive ' //&
@@ -9073,6 +9121,6 @@
 
     USE netcdf_data_input_mod,                                                                     &
-        ONLY:  check_existence, chem_emis_att_type, chem_emis_val_type, get_attribute,             &
-               get_variable, inquire_num_variables, inquire_variable_names,                        &
+        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
+               inquire_num_variables, inquire_variable_names,                                      &
                netcdf_data_input_get_dimension_length, open_read_file
 
@@ -9098,7 +9146,4 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_factor  !< emission time factor
 
-    TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
-    TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
-
 !
 !-- Reset surface fluxes
@@ -9108,15 +9153,16 @@
 
 #if defined( __netcdf )
+!
+!-- Check existence of PIDS_CHEM file
+    INQUIRE( FILE = 'PIDS_CHEM' // TRIM( coupling_char ), EXIST = netcdf_extend )
+    IF ( .NOT. netcdf_extend )  THEN
+       message_string = 'Input file PIDS_CHEM' //  TRIM( coupling_char ) // ' missing!'
+       CALL message( 'salsa_gas_emission_setup', 'PA0640', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Open file in read-only mode
+    CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_chem )
+
     IF ( init )  THEN
-!
-!--    Check existence of PIDS_CHEM file
-       INQUIRE( FILE = 'PIDS_CHEM' // TRIM( coupling_char ), EXIST = netcdf_extend )
-       IF ( .NOT. netcdf_extend )  THEN
-          message_string = 'Input file PIDS_CHEM' //  TRIM( coupling_char ) // ' missing!'
-          CALL message( 'salsa_gas_emission_setup', 'PA0640', 1, 2, 0, 6, 0 )
-       ENDIF
-!
-!--    Open file in read-only mode
-       CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_chem )
 !
 !--    Read the index and name of chemical components
@@ -9301,4 +9347,8 @@
        DEALLOCATE( chem_emission%default_emission_data )
     ENDIF
+!
+!-- Close input file
+    CALL close_input_file( id_chem )
+
 #else
     message_string = 'salsa_emission_mode = "read_from_file", but preprocessor directive ' //   &
@@ -9347,4 +9397,5 @@
        TYPE(surf_type), INTENT(inout) :: surface  !< respective surface type
 
+       conv = 1.0_wp
        use_time_fac = PRESENT( time_fac )
 
@@ -9765,5 +9816,5 @@
              ENDDO
 
-          CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
              IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
                 found_index = get_index( prtcl, TRIM( variable(3:) ) )
@@ -9786,4 +9837,17 @@
                 ENDDO
              ENDIF
+
+          CASE ( 's_H2O' )
+             found_index = get_index( prtcl,'H2O' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
+                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                           &
+                                                 aerosol_mass(ic)%conc(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
 
           CASE DEFAULT
@@ -9885,7 +9949,6 @@
              ENDDO
 
-          CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
              IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
-                found_index = get_index( prtcl, TRIM( variable(3:) ) )
                 IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_bc_av
                 IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_du_av
@@ -9904,4 +9967,15 @@
                 ENDDO
              ENDIF
+
+          CASE ( 's_H2O' )
+             to_be_resorted => s_h2o_av
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                              &
+                                              REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
 
        END SELECT
@@ -10212,5 +10286,5 @@
           IF ( mode == 'xy' )  grid = 'zu'
 
-       CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+       CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
           vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) )
           IF ( is_used( prtcl, vari ) )  THEN
@@ -10252,4 +10326,33 @@
           IF ( mode == 'xy' )  grid = 'zu'
 
+       CASE ( 's_H2O' )
+          found_index = get_index( prtcl, 'H2O' )
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      temp_bin = 0.0_wp
+                      DO  ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol
+                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
+                      ENDDO
+                      local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             to_be_resorted => s_h2o_av
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
+                                               KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( mode == 'xy' )  grid = 'zu'
+
        CASE DEFAULT
           found = .FALSE.
@@ -10534,5 +10637,5 @@
           ENDIF
 
-       CASE ( 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+       CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
           IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
              found_index = get_index( prtcl, TRIM( variable(3:) ) )
@@ -10567,4 +10670,31 @@
                 ENDDO
              ENDIF
+          ENDIF
+
+       CASE ( 's_H2O' )
+          found_index = get_index( prtcl, 'H2O' )
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      temp_bin = 0.0_wp
+                      DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
+                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
+                      ENDDO
+                      local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             to_be_resorted => s_h2o_av
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,            &
+                                               KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
+                   ENDDO
+                ENDDO
+             ENDDO
           ENDIF
 
@@ -10962,4 +11092,43 @@
           ENDIF
 
+       CASE ( 's_H2O' )
+          IF ( av == 0 )  THEN
+             found_index = get_index( prtcl, 'H2O' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nz_do3d
+                      temp_bin = 0.0_wp
+                      DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
+                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
+                      ENDDO
+                      tend(k,j,i) = temp_bin
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. mask_surface(mid) )  THEN
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i),    &
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),         &
+                                                 mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+          ELSE
+             to_be_resorted => s_h2o_av
+          ENDIF
+
        CASE DEFAULT
           found = .FALSE.
@@ -10999,3 +11168,98 @@
  END SUBROUTINE salsa_data_output_mask
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Creates index tables for different (aerosol) components 
+!------------------------------------------------------------------------------!
+ SUBROUTINE component_index_constructor( self, ncomp, nlist, listcomp )
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  ii  !<
+    INTEGER(iwp) ::  jj  !<
+
+    INTEGER(iwp), INTENT(in) ::  nlist ! < Maximum number of components
+
+    INTEGER(iwp), INTENT(inout) ::  ncomp  !< Number of components
+
+    CHARACTER(len=3), INTENT(in) ::  listcomp(nlist)  !< List cof component names
+
+    TYPE(component_index), INTENT(inout) ::  self  !< Object containing the indices of different
+                                                   !< aerosol components
+
+    ncomp = 0
+
+    DO WHILE ( listcomp(ncomp+1) /= '  ' .AND. ncomp < nlist )
+       ncomp = ncomp + 1
+    ENDDO
+
+    self%ncomp = ncomp
+    ALLOCATE( self%ind(ncomp), self%comp(ncomp) )
+
+    DO  ii = 1, ncomp
+       self%ind(ii) = ii
+    ENDDO
+
+    jj = 1
+    DO  ii = 1, nlist
+       IF ( listcomp(ii) == '') CYCLE
+       self%comp(jj) = listcomp(ii)
+       jj = jj + 1
+    ENDDO
+
+ END SUBROUTINE component_index_constructor
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Gives the index of a component in the component list
+!------------------------------------------------------------------------------!
+ INTEGER FUNCTION get_index( self, incomp )
+
+    IMPLICIT NONE
+
+    CHARACTER(len=*), INTENT(in) ::  incomp !< Component name
+
+    INTEGER(iwp) ::  ii  !< index
+
+    TYPE(component_index), INTENT(in) ::  self  !< Object containing the indices of different
+                                                !< aerosol components
+    IF ( ANY( self%comp == incomp ) )  THEN
+       ii = 1
+       DO WHILE ( (self%comp(ii) /= incomp) )
+          ii = ii + 1
+       ENDDO
+       get_index = ii
+    ELSEIF ( incomp == 'H2O' )  THEN
+       get_index = self%ncomp + 1
+    ELSE
+       WRITE( message_string, * ) 'Incorrect component name given!'
+       CALL message( 'get_index', 'PA0591', 1, 2, 0, 6, 0 )
+    ENDIF
+
+ END FUNCTION get_index
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Tells if the (aerosol) component is being used in the simulation
+!------------------------------------------------------------------------------!
+ LOGICAL FUNCTION is_used( self, icomp )
+
+    IMPLICIT NONE
+
+    CHARACTER(len=*), INTENT(in) ::  icomp !< Component name
+
+    TYPE(component_index), INTENT(in) ::  self  !< Object containing the indices of different
+                                                !< aerosol components
+
+    IF ( ANY(self%comp == icomp) ) THEN
+       is_used = .TRUE.
+    ELSE
+       is_used = .FALSE.
+    ENDIF
+
+ END FUNCTION
+
  END MODULE salsa_mod
