Changeset 3956 for palm

Timestamp:

May 7, 2019 12:32:52 PM (6 years ago)

Author:

monakurppa

Message:

Remove salsa calls from prognostic_equations and correct a bug in the salsa deposition for urban and land surface models

Location:

palm/trunk/SOURCE

Files:

• module_interface.f90 (modified) (8 diffs)
• prognostic_equations.f90 (modified) (6 diffs)
• salsa_mod.f90 (modified) (40 diffs)

palm/trunk/SOURCE/module_interface.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! - Added calls for salsa_non_advective_processes and 
+!   salsa_exchange_horiz_bounds
+! - Moved the call for salsa_data_output_2d/3d before that of
+!   radiation_data_output_2d/3d. radiation_data_output_2d/3d tries to read a
+!   salsa output variable and encounters a segmentation fault for "Ntot" due
+!   to the shortoutput name
+! 
+! 3931 2019-04-24 16:34:28Z schwenkel
 ! Changed non_transport_physics to non_advective_processes
 ! 
@@ -158 +166 @@
                bcm_actions,                                                    &
                bcm_non_advective_processes,                                    &
-               bcm_exchange_horiz,                                             &               
+               bcm_exchange_horiz,                                             &
                bcm_prognostic_equations,                                       &
                bcm_swap_timelevel,                                             &
@@ -303 +311 @@
                salsa_header,                                                   &
                salsa_actions,                                                  &
+               salsa_non_advective_processes,                                  &
+               salsa_exchange_horiz_bounds,                                    &
                salsa_prognostic_equations,                                     &
                salsa_swap_timelevel,                                           &
@@ -963 +973 @@
     IF ( bulk_cloud_model    )  CALL bcm_non_advective_processes()
     IF ( air_chemistry       )  CALL chem_non_advective_processes()
+    IF ( salsa               )  CALL salsa_non_advective_processes()
 
 
@@ -980 +991 @@
 
 
-    IF ( bulk_cloud_model    )  CALL bcm_non_advective_processes( i, j )    
+    IF ( bulk_cloud_model    )  CALL bcm_non_advective_processes( i, j )
     IF ( air_chemistry       )  CALL chem_non_advective_processes( i, j )
+    IF ( salsa               )  CALL salsa_non_advective_processes( i, j )
 
 
@@ -996 +1008 @@
     IF ( bulk_cloud_model    )  CALL bcm_exchange_horiz()
     IF ( air_chemistry       )  CALL chem_exchange_horiz_bounds()
+    IF ( salsa               )  CALL salsa_exchange_horiz_bounds()
 
  END SUBROUTINE module_interface_exchange_horiz
@@ -1149 +1162 @@
     ENDIF
 
+    IF ( .NOT. found  .AND.  salsa )  THEN
+       CALL salsa_data_output_2d(                                              &
+               av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do&
+            )
+    ENDIF
+
     IF ( .NOT. found  .AND.  radiation )  THEN
        CALL radiation_data_output_2d(                                          &
-               av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do&
-            )
-    ENDIF
-
-    IF ( .NOT. found  .AND.  salsa )  THEN
-       CALL salsa_data_output_2d(                                              &
                av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do&
             )
@@ -1226 +1239 @@
     ENDIF
 
+    IF ( .NOT. found  .AND.  salsa )  THEN
+       CALL salsa_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
+       resorted = .TRUE.
+    ENDIF
+
     IF ( .NOT. found  .AND.  radiation )  THEN
        CALL radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
-       resorted = .TRUE.
-    ENDIF
-
-    IF ( .NOT. found  .AND.  salsa )  THEN
-       CALL salsa_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
        resorted = .TRUE.
     ENDIF

palm/trunk/SOURCE/prognostic_equations.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Removed salsa calls.
+! 
+! 3931 2019-04-24 16:34:28Z schwenkel
 ! Correct/complete module_interface introduction for chemistry model
 ! 
@@ -464 +467 @@
 #endif
 
-    USE salsa_mod,                                                             &
-        ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time,        &
-               nbins_aerosol, ncomponents_mass, ngases_salsa,                  &
-               salsa_boundary_conds, salsa_diagnostics, salsa_driver,          &
-               salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
-
     USE statistics,                                                            &
         ONLY:  hom
@@ -518 +515 @@
     INTEGER(iwp) ::  i                   !<
     INTEGER(iwp) ::  i_omp_start         !<
-    INTEGER(iwp) ::  ib                  !< index for aerosol size bins (salsa)
-    INTEGER(iwp) ::  ic                  !< index for chemical compounds (salsa)
-    INTEGER(iwp) ::  icc                 !< additional index for chemical compounds (salsa)
-    INTEGER(iwp) ::  ig                  !< index for gaseous compounds (salsa)
     INTEGER(iwp) ::  j                   !<
     INTEGER(iwp) ::  k                   !<
@@ -550 +543 @@
 !
 !-- Module Inferface for exchange horiz after non_advective_processes but before
-!-- advection. Therefore, non_advective_processes must not run for ghost points. 
+!-- advection. Therefore, non_advective_processes must not run for ghost points.
+    !$OMP END PARALLEL
     CALL module_interface_exchange_horiz()
-    !$OMP END PARALLEL
-
-!
-!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time 
-!-- step. The exchange of ghost points is required after this update of the 
-!-- concentrations of aerosol number and mass 
-    IF ( salsa )  THEN
-
-       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
-          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
-          THEN
-             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
-             !$OMP PARALLEL PRIVATE (i,j)
-!
-!--          Call salsa processes
-             !$OMP DO
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   CALL salsa_diagnostics( i, j )
-                   CALL salsa_driver( i, j, 3 )
-                   CALL salsa_diagnostics( i, j )
-                ENDDO
-             ENDDO
-             !$OMP END PARALLEL
-
-             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
-
-             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
-!
-!--          Exchange ghost points and decycle if needed.
-             DO  ib = 1, nbins_aerosol
-                CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
-                CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
-                                           aerosol_number(ib)%init )
-                DO  ic = 1, ncomponents_mass
-                   icc = ( ic - 1 ) * nbins_aerosol + ib
-                   CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
-                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
-                                              aerosol_mass(icc)%init )
-                ENDDO
-             ENDDO
-
-             IF ( .NOT. salsa_gases_from_chem )  THEN
-                DO  ig = 1, ngases_salsa
-                   CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
-                   CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
-                                              salsa_gas(ig)%init )
-                ENDDO
-             ENDIF
-             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
-             last_salsa_time = time_since_reference_point
-
-          ENDIF
-
-       ENDIF
-
-    ENDIF
-
 !
 !-- Loop over all prognostic equations
@@ -1153 +1089 @@
 
     INTEGER(iwp) ::  i     !<
-    INTEGER(iwp) ::  ib    !< index for aerosol size bins (salsa)
-    INTEGER(iwp) ::  ic    !< index for chemical compounds (salsa)
-    INTEGER(iwp) ::  icc   !< additional index for chemical compounds (salsa)
-    INTEGER(iwp) ::  ig    !< index for gaseous compounds (salsa)
     INTEGER(iwp) ::  j     !<
     INTEGER(iwp) ::  k     !<
@@ -1168 +1100 @@
     ENDIF
 !
-!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time 
-!-- step. The exchange of ghost points is required after this update of the 
-!-- concentrations of aerosol number and mass 
-    IF ( salsa )  THEN
-       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
-          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
-          THEN
-             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
-             !$OMP PARALLEL PRIVATE (i,j)
-             !$OMP DO
-!
-!--          Call salsa processes
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   CALL salsa_diagnostics( i, j )
-                   CALL salsa_driver( i, j, 3 )
-                   CALL salsa_diagnostics( i, j )
-                ENDDO
-             ENDDO
-             !$OMP END PARALLEL
-
-             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
-             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
-!
-!--          Exchange ghost points and decycle if needed.
-             DO  ib = 1, nbins_aerosol
-                CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
-                CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
-                                           aerosol_number(ib)%init )
-                DO  ic = 1, ncomponents_mass
-                   icc = ( ic - 1 ) * nbins_aerosol + ib
-                   CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
-                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
-                                              aerosol_mass(icc)%init )
-                ENDDO
-             ENDDO
-             IF ( .NOT. salsa_gases_from_chem )  THEN
-                DO  ig = 1, ngases_salsa
-                   CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
-                   CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
-                                              salsa_gas(ig)%init )
-                ENDDO
-             ENDIF
-             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
-             last_salsa_time = time_since_reference_point
-          ENDIF
-       ENDIF
-    ENDIF
-
-!
 !-- Calculate non advective processes for all other modules
     CALL module_interface_non_advective_processes
 !
 !-- Module Inferface for exchange horiz after non_advective_processes but before
-!-- advection. Therefore, non_advective_processes must not run for ghost points.     
+!-- advection. Therefore, non_advective_processes must not run for ghost points.
     CALL module_interface_exchange_horiz()
-    
 !
 !-- u-velocity component

palm/trunk/SOURCE/salsa_mod.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! - Conceptual bug in depo_surf correct for urban and land surface model
+! - Subroutine salsa_tendency_ij optimized.
+! - Interfaces salsa_non_advective_processes and salsa_exchange_horiz_bounds
+!   created. These are now called in module_interface.
+!   salsa_exchange_horiz_bounds after calling salsa_driver only when needed
+!   (i.e. every dt_salsa).
+! 
+! 3924 2019-04-23 09:33:06Z monakurppa
 ! Correct a bug introduced by the previous update.
 ! 
@@ -173 +181 @@
 !-- Local constants:
     INTEGER(iwp), PARAMETER ::  luc_urban = 15     !< default landuse type for urban
-    INTEGER(iwp), PARAMETER ::  ngases_salsa   = 5 !< total number of gaseous tracers:
+    INTEGER(iwp), PARAMETER ::  ngases_salsa  = 5  !< total number of gaseous tracers:
                                                    !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV
                                                    !< (non-volatile OC), 5 = OCSV (semi-volatile)
-    INTEGER(iwp), PARAMETER ::  nmod    = 7  !< number of modes for initialising the aerosol size
+    INTEGER(iwp), PARAMETER ::  nmod = 7     !< number of modes for initialising the aerosol size
                                              !< distribution
-    INTEGER(iwp), PARAMETER ::  nreg    = 2  !< Number of main size subranges
+    INTEGER(iwp), PARAMETER ::  nreg = 2     !< Number of main size subranges
     INTEGER(iwp), PARAMETER ::  maxspec = 7  !< Max. number of aerosol species
     INTEGER(iwp), PARAMETER ::  season = 1   !< For dry depostion by Zhang et al.: 1 = summer,
@@ -208 +216 @@
 !
 !-- Molar masses in kg/mol
-    REAL(wp), PARAMETER ::  ambc   = 12.0E-3_wp     !< black carbon (BC)
-    REAL(wp), PARAMETER ::  amdair = 28.970E-3_wp   !< dry air
-    REAL(wp), PARAMETER ::  amdu   = 100.E-3_wp     !< mineral dust
-    REAL(wp), PARAMETER ::  amh2o  = 18.0154E-3_wp  !< H2O
-    REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp  !< H2SO4
-    REAL(wp), PARAMETER ::  amhno3 = 63.01E-3_wp    !< HNO3
-    REAL(wp), PARAMETER ::  amn2o  = 44.013E-3_wp   !< N2O
-    REAL(wp), PARAMETER ::  amnh3  = 17.031E-3_wp   !< NH3
-    REAL(wp), PARAMETER ::  amo2   = 31.9988E-3_wp  !< O2
-    REAL(wp), PARAMETER ::  amo3   = 47.998E-3_wp   !< O3
-    REAL(wp), PARAMETER ::  amoc   = 150.E-3_wp     !< organic carbon (OC)
-    REAL(wp), PARAMETER ::  amss   = 58.44E-3_wp    !< sea salt (NaCl)
+    REAL(wp), PARAMETER ::  ambc     = 12.0E-3_wp     !< black carbon (BC)
+    REAL(wp), PARAMETER ::  amdair   = 28.970E-3_wp   !< dry air
+    REAL(wp), PARAMETER ::  amdu     = 100.E-3_wp     !< mineral dust
+    REAL(wp), PARAMETER ::  amh2o    = 18.0154E-3_wp  !< H2O
+    REAL(wp), PARAMETER ::  amh2so4  = 98.06E-3_wp    !< H2SO4
+    REAL(wp), PARAMETER ::  amhno3   = 63.01E-3_wp    !< HNO3
+    REAL(wp), PARAMETER ::  amn2o    = 44.013E-3_wp   !< N2O
+    REAL(wp), PARAMETER ::  amnh3    = 17.031E-3_wp   !< NH3
+    REAL(wp), PARAMETER ::  amo2     = 31.9988E-3_wp  !< O2
+    REAL(wp), PARAMETER ::  amo3     = 47.998E-3_wp   !< O3
+    REAL(wp), PARAMETER ::  amoc     = 150.E-3_wp     !< organic carbon (OC)
+    REAL(wp), PARAMETER ::  amss     = 58.44E-3_wp    !< sea salt (NaCl)
 !
 !-- Densities in kg/m3
@@ -243 +251 @@
 !-- C_B, C_IN, C_IM, beta_IM and C_IT for each land use category (15, as in Zhang et al. (2001))
     REAL(wp), DIMENSION(1:15), PARAMETER :: l_p10 = &
-        (/0.15, 4.0, 0.15, 3.0, 3.0, 0.5, 3.0, -99., 0.5, 2.0, 1.0, -99., -99., -99., 3.0/)
+        (/0.15,  4.0,   0.15,  3.0,   3.0,   0.5,   3.0,   -99., 0.5,  2.0,  1.0,   -99., -99., -99., 3.0/)
     REAL(wp), DIMENSION(1:15), PARAMETER :: c_b_p10 = &
-        (/0.887, 1.262, 0.887, 1.262, 1.262, 0.996, 0.996, -99., 0.7, 0.93, 0.996, -99., -99., -99., 1.262/)
+        (/0.887, 1.262, 0.887, 1.262, 1.262, 0.996, 0.996, -99., 0.7,  0.93, 0.996, -99., -99., -99., 1.262/)
     REAL(wp), DIMENSION(1:15), PARAMETER :: c_in_p10 = &
-        (/0.81, 0.216, 0.81, 0.216, 0.216, 0.191, 0.162, -99., 0.7, 0.14, 0.162, -99., -99., -99., 0.216/)
+        (/0.81,  0.216, 0.81,  0.216, 0.216, 0.191, 0.162, -99., 0.7,  0.14, 0.162, -99., -99., -99., 0.216/)
     REAL(wp), DIMENSION(1:15), PARAMETER :: c_im_p10 = &
-        (/0.162, 0.13, 0.162, 0.13, 0.13, 0.191, 0.081, -99., 0.191, 0.086, 0.081, -99., -99., -99., 0.13/)
+        (/0.162, 0.13,  0.162, 0.13,  0.13,  0.191, 0.081, -99., 0.191,0.086,0.081, -99., -99., -99., 0.13/)
     REAL(wp), DIMENSION(1:15), PARAMETER :: beta_im_p10 = &
-        (/0.6, 0.47, 0.6, 0.47, 0.47, 0.47, 0.47, -99., 0.6, 0.47, 0.47, -99., -99., -99., 0.47/)
+        (/0.6,   0.47,  0.6,   0.47,  0.47,  0.47,  0.47,  -99., 0.6,  0.47, 0.47,  -99., -99., -99., 0.47/)
     REAL(wp), DIMENSION(1:15), PARAMETER :: c_it_p10 = &
-        (/0.0, 0.056, 0.0, 0.056, 0.056, 0.042, 0.056, -99., 0.042, 0.014, 0.056, -99., -99., -99., 0.056/)
+        (/0.0,   0.056, 0.0,   0.056, 0.056, 0.042, 0.056, -99., 0.042,0.014,0.056, -99., -99., -99., 0.056/)
 !
 !-- Constants for the dry deposition model by Zhang et al. (2001): 
@@ -307 +315 @@
                                    (/'SO4','   ','   ','   ','   ','   ','   '/)
 
+    INTEGER(iwp) ::  depo_pcm_par_num = 1   !< parametrisation type: 1=zhang2001, 2=petroff2010
     INTEGER(iwp) ::  depo_pcm_type_num = 0  !< index for the dry deposition type on the plant canopy
+    INTEGER(iwp) ::  depo_surf_par_num = 1  !< parametrisation type: 1=zhang2001, 2=petroff2010
     INTEGER(iwp) ::  dots_salsa = 0         !< starting index for salsa-timeseries
     INTEGER(iwp) ::  end_subrange_1a = 1    !< last index for bin subrange 1a
@@ -364 +374 @@
 !
 !-- SALSA switches:
-    LOGICAL ::  advect_particle_water = .TRUE.     !< advect water concentration of particles
-    LOGICAL ::  decycle_lr            = .FALSE.    !< Undo cyclic boundary conditions: left and right
-    LOGICAL ::  decycle_ns            = .FALSE.    !< north and south boundaries
-    LOGICAL ::  include_emission      = .FALSE.    !< include or not emissions
-    LOGICAL ::  feedback_to_palm      = .FALSE.    !< allow feedback due to condensation of H2O
-    LOGICAL ::  nest_salsa            = .FALSE.    !< apply nesting for salsa
-    LOGICAL ::  no_insoluble          = .FALSE.    !< Switch to exclude insoluble chemical components
-    LOGICAL ::  read_restart_data_salsa = .FALSE.  !< read restart data for salsa
-    LOGICAL ::  salsa_gases_from_chem = .FALSE.    !< Transfer the gaseous components to SALSA from
-                                                   !< from chemistry model
-    LOGICAL ::  van_der_waals_coagc   = .FALSE.    !< Enhancement of coagulation kernel by van der
+    LOGICAL ::  advect_particle_water   = .TRUE.   !< Advect water concentration of particles
+    LOGICAL ::  decycle_lr              = .FALSE.  !< Undo cyclic boundaries: left and right
+    LOGICAL ::  decycle_ns              = .FALSE.  !< Undo cyclic boundaries: north and south
+    LOGICAL ::  include_emission        = .FALSE.  !< Include or not emissions
+    LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
+    LOGICAL ::  nest_salsa              = .FALSE.  !< Apply nesting for salsa
+    LOGICAL ::  no_insoluble            = .FALSE.  !< Exclude insoluble chemical components
+    LOGICAL ::  read_restart_data_salsa = .FALSE.  !< Read restart data for salsa
+    LOGICAL ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA from
+                                                   !< the chemistry model
+    LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Enhancement of coagulation kernel by van der
                                                    !< Waals and viscous forces
-    LOGICAL ::  write_binary_salsa    = .FALSE.    !< read binary for salsa
+    LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
 !
 !-- Process switches: nl* is read from the NAMELIST and is NOT changed.
@@ -455 +465 @@
 !-- SALSA derived datatypes:
 !
-!-- For matching LSM and the deposition module surface types
-    TYPE match_lsm_depo
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match
-    END TYPE match_lsm_depo
+!-- For matching LSM and USM surface types and the deposition module surface types
+    TYPE match_surface
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_lupg  !< index for pavement / green roofs
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luvw  !< index for vegetation / walls
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_luww  !< index for water / windows
+    END TYPE match_surface
 !
 !-- Aerosol emission data attributes
@@ -598 +610 @@
     TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
 
-    TYPE(match_lsm_depo) ::  lsm_to_depo_h
-
-    TYPE(match_lsm_depo), DIMENSION(0:3) ::  lsm_to_depo_v
+    TYPE(match_surface) ::  lsm_to_depo_h  !< to match the deposition module and horizontal LSM surfaces
+    TYPE(match_surface) ::  usm_to_depo_h  !< to match the deposition module and horizontal USM surfaces
+
+    TYPE(match_surface), DIMENSION(0:3) ::  lsm_to_depo_v  !< to match the deposition mod. and vertical LSM surfaces
+    TYPE(match_surface), DIMENSION(0:3) ::  usm_to_depo_v  !< to match the deposition mod. and vertical USM surfaces
 !
 !-- SALSA variables: as x = x(k,j,i,bin). 
@@ -752 +766 @@
     END INTERFACE salsa_actions
 !
+!-- Non-advective processes (i.e. aerosol dynamic reactions) for salsa variables
+    INTERFACE salsa_non_advective_processes
+       MODULE PROCEDURE salsa_non_advective_processes
+       MODULE PROCEDURE salsa_non_advective_processes_ij
+    END INTERFACE salsa_non_advective_processes
+!
+!-- Exchange horiz for salsa variables
+    INTERFACE salsa_exchange_horiz_bounds
+       MODULE PROCEDURE salsa_exchange_horiz_bounds
+    END INTERFACE salsa_exchange_horiz_bounds
+!
 !-- Prognostics equations for salsa variables
     INTERFACE salsa_prognostic_equations
@@ -775 +800 @@
            salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin,        &
            salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local,     &
-           salsa_actions
+           salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds
 !
 !-- Public parameters, constants and initial values
@@ -2162 +2187 @@
 
     USE surface_mod,                                                                               &
-        ONLY:  surf_lsm_h, surf_lsm_v
+        ONLY:  surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
     IMPLICIT NONE
@@ -2168 +2193 @@
     INTEGER(iwp) ::  l  !< loop index for vertical surfaces
 
+    LOGICAL :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM surfaces)
+
+    IF ( depo_pcm_par == 'zhang2001' )  THEN
+       depo_pcm_par_num = 1
+    ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
+       depo_pcm_par_num = 2
+    ENDIF
+
+    IF ( depo_surf_par == 'zhang2001' )  THEN
+       depo_surf_par_num = 1
+    ELSEIF ( depo_surf_par == 'petroff2010' )  THEN
+       depo_surf_par_num = 2
+    ENDIF
+!
+!-- LSM: Pavement, vegetation and water
     IF ( nldepo_surf  .AND.  land_surface )  THEN
-
-       ALLOCATE( lsm_to_depo_h%match(1:surf_lsm_h%ns) )
-       lsm_to_depo_h%match = 0
-       CALL match_lsm_zhang( surf_lsm_h, lsm_to_depo_h%match )
-
+       match_lsm = .TRUE.
+       ALLOCATE( lsm_to_depo_h%match_lupg(1:surf_lsm_h%ns),                                         &
+                 lsm_to_depo_h%match_luvw(1:surf_lsm_h%ns),                                         &
+                 lsm_to_depo_h%match_luww(1:surf_lsm_h%ns) )
+       lsm_to_depo_h%match_lupg = 0
+       lsm_to_depo_h%match_luvw = 0
+       lsm_to_depo_h%match_luww = 0
+       CALL match_sm_zhang( surf_lsm_h, lsm_to_depo_h%match_lupg, lsm_to_depo_h%match_luvw,        &
+                            lsm_to_depo_h%match_luww, match_lsm )
        DO  l = 0, 3
-          ALLOCATE( lsm_to_depo_v(l)%match(1:surf_lsm_v(l)%ns) )
-          lsm_to_depo_v(l)%match = 0
-          CALL match_lsm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match )
+          ALLOCATE( lsm_to_depo_v(l)%match_lupg(1:surf_lsm_v(l)%ns),                               &
+                    lsm_to_depo_v(l)%match_luvw(1:surf_lsm_v(l)%ns),                               &
+                    lsm_to_depo_v(l)%match_luww(1:surf_lsm_v(l)%ns) )
+          lsm_to_depo_v(l)%match_lupg = 0
+          lsm_to_depo_v(l)%match_luvw = 0
+          lsm_to_depo_v(l)%match_luww = 0
+          CALL match_sm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match_lupg,                         &
+                               lsm_to_depo_v(l)%match_luvw, lsm_to_depo_v(l)%match_luww, match_lsm )
+       ENDDO
+    ENDIF
+!
+!-- USM: Green roofs/walls, wall surfaces and windows
+    IF ( nldepo_surf  .AND.  urban_surface )  THEN
+       match_lsm = .FALSE.
+       ALLOCATE( usm_to_depo_h%match_lupg(1:surf_usm_h%ns),                                        &
+                 usm_to_depo_h%match_luvw(1:surf_usm_h%ns),                                        &
+                 usm_to_depo_h%match_luww(1:surf_usm_h%ns) )
+       usm_to_depo_h%match_lupg = 0
+       usm_to_depo_h%match_luvw = 0
+       usm_to_depo_h%match_luww = 0
+       CALL match_sm_zhang( surf_usm_h, usm_to_depo_h%match_lupg, usm_to_depo_h%match_luvw,        &
+                            usm_to_depo_h%match_luww, match_lsm )
+       DO  l = 0, 3
+          ALLOCATE( usm_to_depo_v(l)%match_lupg(1:surf_usm_v(l)%ns),                               &
+                    usm_to_depo_v(l)%match_luvw(1:surf_usm_v(l)%ns),                               &
+                    usm_to_depo_v(l)%match_luww(1:surf_usm_v(l)%ns) )
+          usm_to_depo_v(l)%match_lupg = 0
+          usm_to_depo_v(l)%match_luvw = 0
+          usm_to_depo_v(l)%match_luww = 0
+          CALL match_sm_zhang( surf_usm_v(l), usm_to_depo_v(l)%match_lupg,                         &
+                               usm_to_depo_v(l)%match_luvw, usm_to_depo_v(l)%match_luww, match_lsm )
        ENDDO
     ENDIF
@@ -2204 +2276 @@
 !> Match the surface types in PALM and Zhang et al. 2001 deposition module
 !------------------------------------------------------------------------------!
- SUBROUTINE match_lsm_zhang( surf, match_array )
+ SUBROUTINE match_sm_zhang( surf, match_pav_green, match_veg_wall, match_wat_win, match_lsm )
 
     USE surface_mod,                                                           &
@@ -2211 +2283 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  m                !< index for surface elements
-    INTEGER(iwp) ::  pav_type_palm    !< pavement type in PALM
-    INTEGER(iwp) ::  vege_type_palm   !< vegetation type in PALM
-    INTEGER(iwp) ::  water_type_palm  !< water type in PALM
-
-    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_array !< array matching
-                                                              !< the surface types
+    INTEGER(iwp) ::  m              !< index for surface elements
+    INTEGER(iwp) ::  pav_type_palm  !< pavement / green wall type in PALM
+    INTEGER(iwp) ::  veg_type_palm  !< vegetation / wall type in PALM
+    INTEGER(iwp) ::  wat_type_palm  !< water / window type in PALM
+
+    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_pav_green  !<  matching pavement/green walls
+    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_veg_wall   !<  matching vegetation/walls
+    INTEGER(iwp), DIMENSION(:), INTENT(inout) ::  match_wat_win    !<  matching water/windows
+
+    LOGICAL, INTENT(in) :: match_lsm  !< flag to initilise LSM surfaces (if false, initialise USM)
+
     TYPE(surf_type), INTENT(in) :: surf  !< respective surface type
 
     DO  m = 1, surf%ns
-
-       IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
-          vege_type_palm = surf%vegetation_type(m)
-          SELECT CASE ( vege_type_palm )
-             CASE ( 0 )
-                message_string = 'No vegetation type defined.'
-                CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
-             CASE ( 1 )  ! bare soil
-                match_array(m) = 6  ! grass in Z01
-             CASE ( 2 )  ! crops, mixed farming
-                match_array(m) = 7  !  crops, mixed farming Z01
-             CASE ( 3 )  ! short grass
-                match_array(m) = 6  ! grass in Z01
-             CASE ( 4 )  ! evergreen needleleaf trees
-                 match_array(m) = 1  ! evergreen needleleaf trees in Z01
-             CASE ( 5 )  ! deciduous needleleaf trees
-                match_array(m) = 3  ! deciduous needleleaf trees in Z01
-             CASE ( 6 )  ! evergreen broadleaf trees
-                match_array(m) = 2  ! evergreen broadleaf trees in Z01
-             CASE ( 7 )  ! deciduous broadleaf trees
-                match_array(m) = 4  ! deciduous broadleaf trees in Z01
-             CASE ( 8 )  ! tall grass
-                match_array(m) = 6  ! grass in Z01
-             CASE ( 9 )  ! desert
-                match_array(m) = 8  ! desert in Z01
-             CASE ( 10 )  ! tundra
-                match_array(m) = 9  ! tundra in Z01
-             CASE ( 11 )  ! irrigated crops
-                match_array(m) = 7  !  crops, mixed farming Z01
-             CASE ( 12 )  ! semidesert
-                match_array(m) = 8  ! desert in Z01
-             CASE ( 13 )  ! ice caps and glaciers
-                match_array(m) = 12  ! ice cap and glacier in Z01
-             CASE ( 14 )  ! bogs and marshes
-                match_array(m) = 11  ! wetland with plants in Z01
-             CASE ( 15 )  ! evergreen shrubs
-                match_array(m) = 10  ! shrubs and interrupted woodlands in Z01
-             CASE ( 16 )  ! deciduous shrubs
-                match_array(m) = 10  ! shrubs and interrupted woodlands in Z01
-             CASE ( 17 )  ! mixed forest/woodland
-                match_array(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
-             CASE ( 18 )  ! interrupted forest
-                match_array(m) = 10  ! shrubs and interrupted woodlands in Z01
-          END SELECT
-       ENDIF
-
-       IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
-          pav_type_palm = surf%pavement_type(m)
-          IF ( pav_type_palm == 0 )  THEN  ! error
-             message_string = 'No pavement type defined.'
-             CALL message( 'salsa_mod: match_lsm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
-          ELSEIF ( pav_type_palm > 0  .AND.  pav_type_palm <= 15 )  THEN
-             match_array(m) = 15  ! urban in Z01
+       IF ( match_lsm )  THEN
+!
+!--       Vegetation (LSM):
+          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
+             veg_type_palm = surf%vegetation_type(m)
+             SELECT CASE ( veg_type_palm )
+                CASE ( 0 )
+                   message_string = 'No vegetation type defined.'
+                   CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 )
+                CASE ( 1 )  ! bare soil
+                   match_veg_wall(m) = 6  ! grass in Z01
+                CASE ( 2 )  ! crops, mixed farming
+                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
+                CASE ( 3 )  ! short grass
+                   match_veg_wall(m) = 6  ! grass in Z01
+                CASE ( 4 )  ! evergreen needleleaf trees
+                    match_veg_wall(m) = 1  ! evergreen needleleaf trees in Z01
+                CASE ( 5 )  ! deciduous needleleaf trees
+                   match_veg_wall(m) = 3  ! deciduous needleleaf trees in Z01
+                CASE ( 6 )  ! evergreen broadleaf trees
+                   match_veg_wall(m) = 2  ! evergreen broadleaf trees in Z01
+                CASE ( 7 )  ! deciduous broadleaf trees
+                   match_veg_wall(m) = 4  ! deciduous broadleaf trees in Z01
+                CASE ( 8 )  ! tall grass
+                   match_veg_wall(m) = 6  ! grass in Z01
+                CASE ( 9 )  ! desert
+                   match_veg_wall(m) = 8  ! desert in Z01
+                CASE ( 10 )  ! tundra
+                   match_veg_wall(m) = 9  ! tundra in Z01
+                CASE ( 11 )  ! irrigated crops
+                   match_veg_wall(m) = 7  !  crops, mixed farming Z01
+                CASE ( 12 )  ! semidesert
+                   match_veg_wall(m) = 8  ! desert in Z01
+                CASE ( 13 )  ! ice caps and glaciers
+                   match_veg_wall(m) = 12  ! ice cap and glacier in Z01
+                CASE ( 14 )  ! bogs and marshes
+                   match_veg_wall(m) = 11  ! wetland with plants in Z01
+                CASE ( 15 )  ! evergreen shrubs
+                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
+                CASE ( 16 )  ! deciduous shrubs
+                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
+                CASE ( 17 )  ! mixed forest/woodland
+                   match_veg_wall(m) = 5  ! mixed broadleaf and needleleaf trees in Z01
+                CASE ( 18 )  ! interrupted forest
+                   match_veg_wall(m) = 10  ! shrubs and interrupted woodlands in Z01
+             END SELECT
           ENDIF
-       ENDIF
-
-       IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
-          water_type_palm = surf%water_type(m)
-          IF ( water_type_palm == 0 )  THEN  ! error
-             message_string = 'No water type defined.'
-             CALL message( 'salsa_mod: match_lsm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
-          ELSEIF ( water_type_palm == 3 )  THEN
-             match_array(m) = 14  ! ocean in Z01
-          ELSEIF ( water_type_palm == 1  .OR.  water_type_palm == 2 .OR.  water_type_palm == 4     &
-                   .OR.  water_type_palm == 5  )  THEN
-             match_array(m) = 13  ! inland water in Z01
+!
+!--       Pavement (LSM):
+          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
+             pav_type_palm = surf%pavement_type(m)
+             IF ( pav_type_palm == 0 )  THEN  ! error
+                message_string = 'No pavement type defined.'
+                CALL message( 'salsa_mod: match_sm_zhang', 'PA0615', 1, 2, 0, 6, 0 )
+             ELSE
+                match_pav_green(m) = 15  ! urban in Z01
+             ENDIF
           ENDIF
+!
+!--       Water (LSM):
+          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
+             wat_type_palm = surf%water_type(m)
+             IF ( wat_type_palm == 0 )  THEN  ! error
+                message_string = 'No water type defined.'
+                CALL message( 'salsa_mod: match_sm_zhang', 'PA0616', 1, 2, 0, 6, 0 )
+             ELSEIF ( wat_type_palm == 3 )  THEN
+                match_wat_win(m) = 14  ! ocean in Z01
+             ELSEIF ( wat_type_palm == 1  .OR.  wat_type_palm == 2 .OR.  wat_type_palm == 4        &
+                      .OR.  wat_type_palm == 5  )  THEN
+                match_wat_win(m) = 13  ! inland water in Z01
+             ENDIF
+          ENDIF
+       ELSE
+!
+!--       Wall surfaces (USM):
+          IF ( surf%frac(ind_veg_wall,m) > 0 )  THEN
+             match_veg_wall(m) = 15  ! urban in Z01
+          ENDIF
+!
+!--       Green walls and roofs (USM):
+          IF ( surf%frac(ind_pav_green,m) > 0 )  THEN
+             match_pav_green(m) =  6 ! (short) grass in Z01
+          ENDIF
+!
+!--       Windows (USM):
+          IF ( surf%frac(ind_wat_win,m) > 0 )  THEN
+             match_wat_win(m) = 15  ! urban in Z01
+          ENDIF
        ENDIF
 
     ENDDO
 
- END SUBROUTINE match_lsm_zhang
+ END SUBROUTINE match_sm_zhang
 
 !------------------------------------------------------------------------------!
@@ -2931 +3028 @@
        ENDIF
 !
-!--    Tendency of water vapour mixing ratio is obtained from the 
-!--    change in RH during SALSA run. This releases heat and changes pt.
-!--    Assumes no temperature change during SALSA run.
+!--    Tendency of water vapour mixing ratio is obtained from the change in RH during SALSA run.
+!--    This releases heat and changes pt. Assumes no temperature change during SALSA run.
 !--    q = r / (1+r), Euler method for integration
 !
@@ -2951 +3047 @@
           CALL depo_surf( i, j, surf_def_h(0), vd, schmidt_num, kvis, in_u, .TRUE. )
           DO  l = 0, 3
-             CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l )
+             CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE. )
           ENDDO
        ELSE
-          CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE. )
+          CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE., usm_to_depo_h )
           DO  l = 0, 3
-             CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l )
+             CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
+                             usm_to_depo_v(l) )
           ENDDO
-          CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE. )
+          CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE., lsm_to_depo_h )
           DO  l = 0, 3
-             CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l )
+             CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE.,            &
+                             lsm_to_depo_v(l) )
           ENDDO
        ENDIF
@@ -3342 +3440 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  ib     !< loop index
-
-    REAL(wp) ::  avis       !< molecular viscocity of air (kg/(m*s))
-    REAL(wp) ::  Cc         !< Cunningham slip-flow correction factor
-    REAL(wp) ::  Kn         !< Knudsen number
-    REAL(wp) ::  lambda     !< molecular mean free path (m)
-    REAL(wp) ::  mdiff      !< particle diffusivity coefficient
-    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)
+    INTEGER(iwp) ::  ib   !< loop index
+    INTEGER(iwp) ::  ic   !< loop index
+
+    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
+    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
+    REAL(wp) ::  c_interception    !< coefficient for interception
+    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
+    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
+    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors,
+                                   !< Table 3 in Z01
+    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
+    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), INTENT(in) ::  adn    !< air density (kg/m3)
@@ -3368 +3478 @@
 !
 !-- Initialise
-    pdn           = 1500.0_wp    ! default value
-    ustar         = 0.0_wp
+    depo  = 0.0_wp
+    pdn   = 1500.0_wp    ! default value
+    ustar = 0.0_wp
 !
 !-- Molecular viscosity of air (Eq. 4.54)
@@ -3382 +3493 @@
 !-- Mean free path (m) (Eq. 15.24)
     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
@@ -3392 +3516 @@
 !
 !--    Cunningham slip-flow correction (Eq. 15.30)
-       Cc = 1.0_wp + Kn * ( 1.249_wp + 0.42_wp * EXP( -0.87_wp / Kn ) )
+       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 * Cc ) / ( 3.0_wp * pi * avis * zdwet )
+       mdiff = ( abo * tk * beta ) / ( 3.0_wp * pi * avis * zdwet )
 !
 !--    Particle Schmidt number (Eq. 15.36)
@@ -3401 +3525 @@
 !
 !--    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, Cc) )
+       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
           ustar = SQRT( cdc ) * mag_u
-          CALL depo_pcm( paero, ib, vc(ib), mag_u, ustar, kvis, schmidt_num(ib), lad )
+          SELECT CASE ( depo_pcm_par_num )
+
+             CASE ( 1 )   ! Zhang et al. (2001)
+                CALL depo_vel_Z01( vc(ib), ustar, schmidt_num(ib), paero(ib)%dwet, alpha,  gamma,  &
+                                   par_a, depo )
+             CASE ( 2 )   ! Petroff & Zhang (2010)
+                CALL depo_vel_P10( vc(ib), mag_u, ustar, kvis, schmidt_num(ib), paero(ib)%dwet,    &
+                                   par_l, c_brownian_diff, c_interception, c_impaction, beta_im,   &
+                                   c_turb_impaction, depo )
+          END SELECT
+!
+!--       Calculate the change in concentrations
+          paero(ib)%numc = paero(ib)%numc - depo * lad * paero(ib)%numc * dt_salsa
+          DO  ic = 1, maxspec+1
+             paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * lad * paero(ib)%volc(ic) * dt_salsa
+          ENDDO
        ENDIF
     ENDDO
@@ -3415 +3554 @@
 ! Description:
 ! ------------
-!> Calculate change in number and volume concentrations due to deposition on
-!> plant canopy.
-!------------------------------------------------------------------------------!
- SUBROUTINE depo_pcm( paero, ib, vc, mag_u, ustar, kvis_a, schmidt_num, lad )
+!> Calculate deposition velocity (m/s) based on Zhan et al. (2001, case 1).
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE depo_vel_Z01( vc, ustar, schmidt_num, diameter, alpha, gamma, par_a, depo )
 
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  ic      !< loop index
-
-    INTEGER(iwp), INTENT(in) ::  ib  !< loop index
-
-    REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
-    REAL(wp) ::  beta_im           !< parameter for turbulent impaction
-    REAL(wp) ::  depo              !< deposition efficiency
-    REAL(wp) ::  c_brownian_diff   !< coefficient for Brownian diffusion
-    REAL(wp) ::  c_impaction       !< coefficient for inertial impaction
-    REAL(wp) ::  c_interception    !< coefficient for interception
-    REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
-    REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
-    REAL(wp) ::  par_a             !< parameter A for the characteristic radius of collectors, 
-                                   !< Table 3 in Z01
-    REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
     REAL(wp) ::  rs                !< overall quasi-laminar resistance for particles
+    REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
+
+    REAL(wp), INTENT(in) ::  alpha        !< parameter, Table 3 in Z01
+    REAL(wp), INTENT(in) ::  gamma        !< parameter, Table 3 in Z01
+    REAL(wp), INTENT(in) ::  par_a        !< parameter A for the characteristic diameter of
+                                          !< collectors, Table 3 in Z01
+    REAL(wp), INTENT(in) ::  diameter     !< particle diameter
+    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
+    REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
+    REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
+
+    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
+
+    IF ( par_a > 0.0_wp )  THEN
+!
+!--    Initialise
+       rs = 0.0_wp
+!
+!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
+       stokes_num = vc * ustar / ( g * par_a )
+!
+!--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)
+       rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -stokes_num**0.5_wp ) *                &
+                 ( schmidt_num**( -gamma ) + ( stokes_num / ( alpha + stokes_num ) )**2 +          &
+                 0.5_wp * ( diameter / par_a )**2 ) ) )
+
+       depo = rs + vc
+
+    ELSE
+       depo = 0.0_wp
+    ENDIF
+
+ END SUBROUTINE depo_vel_Z01
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate deposition velocity (m/s) based on Petroff & Zhang (2010, case 2).
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE depo_vel_P10( vc, mag_u, ustar, kvis_a, schmidt_num, diameter, par_l, c_brownian_diff, &
+                          c_interception, c_impaction, beta_im, c_turb_impaction, depo )
+
+    IMPLICIT NONE
+
     REAL(wp) ::  stokes_num        !< Stokes number for smooth or bluff surfaces
     REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
@@ -3445 +3614 @@
     REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
 
+    REAL(wp), INTENT(in) ::  beta_im           !< parameter for turbulent impaction
+    REAL(wp), INTENT(in) ::  c_brownian_diff   !< coefficient for Brownian diffusion
+    REAL(wp), INTENT(in) ::  c_impaction       !< coefficient for inertial impaction
+    REAL(wp), INTENT(in) ::  c_interception    !< coefficient for interception
+    REAL(wp), INTENT(in) ::  c_turb_impaction  !< coefficient for turbulent impaction
     REAL(wp), INTENT(in) ::  kvis_a       !< kinematic viscosity of air (m2/s)
-    REAL(wp), INTENT(in) ::  lad          !< leaf area density (m2/m3)
     REAL(wp), INTENT(in) ::  mag_u        !< wind velocity (m/s)
-    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number 
+    REAL(wp), INTENT(in) ::  par_l        !< obstacle characteristic dimension in P10
+    REAL(wp), INTENT(in) ::  diameter       !< particle diameter
+    REAL(wp), INTENT(in) ::  schmidt_num  !< particle Schmidt number
     REAL(wp), INTENT(in) ::  ustar        !< friction velocity (m/s)
     REAL(wp), INTENT(in) ::  vc           !< terminal velocity (m/s)
 
-    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< aerosol properties
-!
-!-- Initialise
-    rs       = 0.0_wp
-    tau_plus = 0.0_wp
-    v_bd     = 0.0_wp
-    v_im     = 0.0_wp
-    v_in     = 0.0_wp
-    v_it     = 0.0_wp
-
-    IF ( depo_pcm_par == 'zhang2001' )  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
-!
-!--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
-       stokes_num = vc * ustar / ( g * par_a )
-!
-!--    The overall quasi-laminar resistance for particles (Zhang et al., Eq. 5)
-       rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * ustar * EXP( -stokes_num**0.5_wp ) *                &
-                 ( schmidt_num**( -gamma ) + ( stokes_num / ( alpha + stokes_num ) )**2 +          &
-                 0.5_wp * ( paero(ib)%dwet / par_a )**2 ) ) )
-
-       depo = ( rs + vc ) * lad
-
-    ELSEIF ( depo_pcm_par == 'petroff2010' )  THEN
-!
-!--    vd = v_BD + v_IN + v_IM + v_IT + vc
-!--    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)
+    REAL(wp), INTENT(inout)  ::  depo     !< deposition efficiency (m/s)
+
+    IF ( par_l > 0.0_wp )  THEN
+!
+!--    Initialise
+       tau_plus = 0.0_wp
+       v_bd     = 0.0_wp
+       v_im     = 0.0_wp
+       v_in     = 0.0_wp
+       v_it     = 0.0_wp
 !
 !--    Stokes number for vegetated surfaces (Seinfeld & Pandis (2006): Eq.19.24)
@@ -3497 +3645 @@
 !
 !--    Brownian diffusion
-       v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                          &
+       v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                             &
               ( mag_u * par_l / kvis_a )**( -0.5_wp )
 !
 !--    Interception
-       v_in = mag_u * c_interception * paero(ib)%dwet / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l /    &
-              paero(ib)%dwet ) )
+       v_in = mag_u * c_interception * diameter / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) )
 !
 !--    Impaction: Petroff (2009) Eq. 18
@@ -3514 +3661 @@
        ENDIF
 
-       depo = ( v_bd + v_in + v_im + v_it + vc ) * lad
-
+       depo = ( v_bd + v_in + v_im + v_it + vc )
+
+    ELSE
+       depo = 0.0_wp
     ENDIF
-!
-!-- Calculate the change in concentrations
-    paero(ib)%numc = paero(ib)%numc - depo * paero(ib)%numc * dt_salsa
-    DO  ic = 1, maxspec+1
-       paero(ib)%volc(ic) = paero(ib)%volc(ic) - depo * paero(ib)%volc(ic) * dt_salsa
-    ENDDO
-
- END SUBROUTINE depo_pcm
+
+ END SUBROUTINE depo_vel_P10
 
 !------------------------------------------------------------------------------!
@@ -3534 +3677 @@
 !        high-resolution simulations)
 !------------------------------------------------------------------------------!
- SUBROUTINE depo_surf( i, j, surf, vc, schmidt_num, kvis, mag_u, norm, l )
-
-    USE arrays_3d,                                                             &
+ SUBROUTINE depo_surf( i, j, surf, vc, schmidt_num, kvis, mag_u, norm, match_array )
+
+    USE arrays_3d,                                                                                 &
         ONLY: rho_air_zw
 
-    USE surface_mod,                                                           &
-        ONLY:  surf_type
+    USE surface_mod,                                                                               &
+        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_type
 
     IMPLICIT NONE
@@ -3552 +3695 @@
     INTEGER(iwp) ::  surf_s  !< Start index of surface elements at (j,i)-gridpoint
 
-    INTEGER(iwp), INTENT(in) ::  i     !< loop index
-    INTEGER(iwp), INTENT(in) ::  j     !< loop index
-
-    INTEGER(iwp), INTENT(in), OPTIONAL ::  l  !< index variable for surface facing
-
-    LOGICAL, INTENT(in) ::  norm      !< to normalise or not
+    INTEGER(iwp), INTENT(in) ::  i  !< loop index
+    INTEGER(iwp), INTENT(in) ::  j  !< loop index
+
+    LOGICAL, INTENT(in) ::  norm   !< to normalise or not
 
     REAL(wp) ::  alpha             !< parameter, Table 3 in Z01
@@ -3565 +3706 @@
     REAL(wp) ::  c_interception    !< coefficient for interception
     REAL(wp) ::  c_turb_impaction  !< coefficient for turbulent impaction
-    REAL(wp) ::  depo              !< deposition efficiency
     REAL(wp) ::  gamma             !< parameter, Table 3 in Z01
     REAL(wp) ::  norm_fac          !< normalisation factor (usually air density)
@@ -3572 +3712 @@
     REAL(wp) ::  par_l             !< obstacle characteristic dimension in P10
     REAL(wp) ::  rs                !< the overall quasi-laminar resistance for particles
-    REAL(wp) ::  stokes_num        !< Stokes number for bluff surface elements
     REAL(wp) ::  tau_plus          !< dimensionless particle relaxation time
     REAL(wp) ::  v_bd              !< deposition velocity due to Brownian diffusion
@@ -3579 +3718 @@
     REAL(wp) ::  v_it              !< deposition velocity due to turbulent impaction
 
+    REAL(wp), DIMENSION(nbins_aerosol) ::  depo      !< deposition efficiency
+    REAL(wp), DIMENSION(nbins_aerosol) ::  depo_sum  !< sum of deposition efficiencies
+
     REAL(wp), DIMENSION(:), INTENT(in) ::  kvis   !< kinematic viscosity of air (m2/s)
     REAL(wp), DIMENSION(:), INTENT(in) ::  mag_u  !< wind velocity (m/s)
@@ -3585 +3727 @@
     REAL(wp), DIMENSION(:,:), INTENT(in) ::  vc            !< terminal velocity (m/s)
 
-    TYPE(surf_type), INTENT(inout) :: surf  !< respective surface type
+    TYPE(match_surface), INTENT(in), OPTIONAL ::  match_array  !< match the deposition module and
+                                                               !< LSM/USM surfaces
+    TYPE(surf_type), INTENT(inout) :: surf                     !< respective surface type
 !
 !-- Initialise
+    depo     = 0.0_wp
+    depo_sum = 0.0_wp
     rs       = 0.0_wp
     surf_s   = surf%start_index(j,i)
@@ -3597 +3743 @@
     v_it     = 0.0_wp
 !
-!-- Model parameters for the land use category. If LSM is applied, import
+!-- Model parameters for the land use category. If LSM or USM is applied, import
 !-- characteristics. Otherwise, apply surface type "urban".
     alpha   = alpha_z01(luc_urban)
@@ -3610 +3756 @@
     c_turb_impaction = c_it_p10(luc_urban)
 
-    DO  m = surf_s, surf_e
-       k = surf%k(m)
-
-       IF ( norm )  THEN
-          norm_fac = rho_air_zw(k)
-          IF ( land_surface )  THEN
-             alpha            = alpha_z01( lsm_to_depo_h%match(m) )
-             beta_im          = beta_im_p10( lsm_to_depo_h%match(m) )
-             c_brownian_diff  = c_b_p10( lsm_to_depo_h%match(m) )
-             c_impaction      = c_im_p10( lsm_to_depo_h%match(m) )
-             c_interception   = c_in_p10( lsm_to_depo_h%match(m) )
-             c_turb_impaction = c_it_p10( lsm_to_depo_h%match(m) )
-             gamma            = gamma_z01( lsm_to_depo_h%match(m) )
-             par_a            = A_z01( lsm_to_depo_h%match(m), season ) * 1.0E-3_wp
-             par_l            = l_p10( lsm_to_depo_h%match(m) ) * 0.01_wp
+
+    IF ( PRESENT( match_array ) )  THEN  ! land or urban surface model
+
+       DO  m = surf_s, surf_e
+
+          k = surf%k(m)
+          norm_fac = 1.0_wp
+
+          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
+
+          IF ( match_array%match_lupg(m) > 0 )  THEN
+             alpha = alpha_z01( match_array%match_lupg(m) )
+             gamma = gamma_z01( match_array%match_lupg(m) )
+             par_a = A_z01( match_array%match_lupg(m), season ) * 1.0E-3_wp
+
+             beta_im          = beta_im_p10( match_array%match_lupg(m) )
+             c_brownian_diff  = c_b_p10( match_array%match_lupg(m) )
+             c_impaction      = c_im_p10( match_array%match_lupg(m) )
+             c_interception   = c_in_p10( match_array%match_lupg(m) )
+             c_turb_impaction = c_it_p10( match_array%match_lupg(m) )
+             par_l            = l_p10( match_array%match_lupg(m) ) * 0.01_wp
+
+             DO  ib = 1, nbins_aerosol
+                IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
+                   CYCLE
+
+                SELECT CASE ( depo_surf_par_num )
+
+                   CASE ( 1 )
+                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
+                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
+                   CASE ( 2 )
+                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
+                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
+                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
+                                         c_turb_impaction, depo(ib) )
+                END SELECT
+             ENDDO
+             depo_sum = depo_sum + surf%frac(ind_pav_green,m) * depo
           ENDIF
-       ELSE
-          norm_fac = 0.0_wp
-          IF ( land_surface )  THEN
-             alpha            = alpha_z01( lsm_to_depo_v(l)%match(m) )
-             beta_im          = beta_im_p10( lsm_to_depo_v(l)%match(m) )
-             c_brownian_diff  = c_b_p10( lsm_to_depo_v(l)%match(m) )
-             c_impaction      = c_im_p10( lsm_to_depo_v(l)%match(m) )
-             c_interception   = c_in_p10( lsm_to_depo_v(l)%match(m) )
-             c_turb_impaction = c_it_p10( lsm_to_depo_v(l)%match(m) )
-             gamma            = gamma_z01( lsm_to_depo_v(l)%match(m) )
-             par_a            = A_z01( lsm_to_depo_v(l)%match(m), season ) * 1.0E-3_wp
-             par_l            = l_p10( lsm_to_depo_v(l)%match(m) ) * 0.01_wp
+
+          IF ( match_array%match_luvw(m) > 0 )  THEN
+             alpha = alpha_z01( match_array%match_luvw(m) )
+             gamma = gamma_z01( match_array%match_luvw(m) )
+             par_a = A_z01( match_array%match_luvw(m), season ) * 1.0E-3_wp
+
+             beta_im          = beta_im_p10( match_array%match_luvw(m) )
+             c_brownian_diff  = c_b_p10( match_array%match_luvw(m) )
+             c_impaction      = c_im_p10( match_array%match_luvw(m) )
+             c_interception   = c_in_p10( match_array%match_luvw(m) )
+             c_turb_impaction = c_it_p10( match_array%match_luvw(m) )
+             par_l            = l_p10( match_array%match_luvw(m) ) * 0.01_wp
+
+             DO  ib = 1, nbins_aerosol
+                IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
+                   CYCLE
+
+                SELECT CASE ( depo_surf_par_num )
+
+                   CASE ( 1 )
+                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
+                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
+                   CASE ( 2 )
+                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
+                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
+                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
+                                         c_turb_impaction, depo(ib) )
+                END SELECT
+             ENDDO
+             depo_sum = depo_sum + surf%frac(ind_veg_wall,m) * depo
           ENDIF
-       ENDIF
-
-       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 ( depo_surf_par == 'zhang2001' )  THEN
-!
-!--          Stokes number for smooth surfaces or surfaces with bluff roughness
-!--          elements (Seinfeld and Pandis, 2nd edition (2006): Eq. 19.23)
-             stokes_num = MAX( 0.01_wp, vc(k+1,ib) * surf%us(m)**2 / ( g * kvis(k+1)  ) )
-!
-!--          The overall quasi-laminar resistance for particles (Eq. 5)
-             rs = MAX( EPSILON( 1.0_wp ), ( 3.0_wp * surf%us(m) * ( schmidt_num(k+1,ib)**( -gamma )&
-                       + ( stokes_num / ( alpha + stokes_num ) )**2 + 0.5_wp * ( ra_dry(k,j,i,ib) /&
-                       par_a )**2 ) * EXP( -stokes_num**0.5_wp ) ) )
-             depo = vc(k+1,ib) + rs
-
-          ELSEIF ( depo_surf_par == 'petroff2010' )  THEN 
-!
-!--          vd = v_BD + v_IN + v_IM + v_IT + vc
-!
-!--          Stokes number for smooth surfaces or surfaces with bluff roughness
-!--          elements (Seinfeld and Pandis, 2nd edition (2006): Eq. 19.23)
-             stokes_num = MAX( 0.01_wp, vc(k+1,ib) * surf%us(m)**2 / ( g *  kvis(k+1) ) )
-!
-!--          Non-dimensional relexation time of the particle on top of canopy
-             tau_plus = vc(k+1,ib) * surf%us(m)**2 / ( kvis(k+1) * g )
-!
-!--          Brownian diffusion
-             v_bd = mag_u(k+1) * c_brownian_diff * schmidt_num(k+1,ib)**( -0.666666_wp ) *         &
-                    ( mag_u(k+1) * par_l / kvis(k+1) )**( -0.5_wp )
-!
-!--          Interception
-             v_in = mag_u(k+1) * c_interception * ra_dry(k,j,i,ib)/ par_l *                        &
-                    ( 2.0_wp + LOG( 2.0_wp * par_l / ra_dry(k,j,i,ib) ) )
-!
-!--          Impaction: Petroff (2009) Eq. 18
-             v_im = mag_u(k+1) * c_impaction * ( stokes_num / ( stokes_num + beta_im ) )**2
-
-             IF ( tau_plus < 20.0_wp )  THEN
-                v_it = 2.5E-3_wp * c_turb_impaction * tau_plus**2
+
+          IF ( match_array%match_luww(m) > 0 )  THEN
+             alpha = alpha_z01( match_array%match_luww(m) )
+             gamma = gamma_z01( match_array%match_luww(m) )
+             par_a = A_z01( match_array%match_luww(m), season ) * 1.0E-3_wp
+
+             beta_im          = beta_im_p10( match_array%match_luww(m) )
+             c_brownian_diff  = c_b_p10( match_array%match_luww(m) )
+             c_impaction      = c_im_p10( match_array%match_luww(m) )
+             c_interception   = c_in_p10( match_array%match_luww(m) )
+             c_turb_impaction = c_it_p10( match_array%match_luww(m) )
+             par_l            = l_p10( match_array%match_luww(m) ) * 0.01_wp
+
+             DO  ib = 1, nbins_aerosol
+                IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp ) &
+                   CYCLE
+
+                SELECT CASE ( depo_surf_par_num )
+
+                   CASE ( 1 )
+                      CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),              &
+                                         ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
+                   CASE ( 2 )
+                      CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),            &
+                                         schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,             &
+                                         c_brownian_diff, c_interception, c_impaction, beta_im,    &
+                                         c_turb_impaction, depo(ib) )
+                END SELECT
+             ENDDO
+             depo_sum = depo_sum + surf%frac(ind_wat_win,m) * depo
+          ENDIF
+
+          DO  ib = 1, nbins_aerosol
+!
+!--          Calculate changes in surface fluxes due to dry deposition
+             IF ( include_emission )  THEN
+                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
+                                   depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
+                                       depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
+                ENDDO  ! ic
              ELSE
-                v_it = c_turb_impaction
+                surf%answs(m,ib) = -depo_sum(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   surf%amsws(m,icc) = -depo_sum(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
+                ENDDO  ! ic
              ENDIF
-             depo =  v_bd + v_in + v_im + v_it + vc(k+1,ib)
-
-          ENDIF
-!
-!--       Calculate changes in surface fluxes due to dry deposition
-          IF ( include_emission )  THEN
-             surf%answs(m,ib) = aerosol_number(ib)%source(j,i) -                                   &
-                                MAX( 0.0_wp, depo * norm_fac * aerosol_number(ib)%conc(k,j,i) )
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) -                                &
-                                    MAX( 0.0_wp, depo *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
-             ENDDO  ! ic
-          ELSE
-             surf%answs(m,ib) = -depo * norm_fac * aerosol_number(ib)%conc(k,j,i)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                surf%amsws(m,icc) = -depo *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
-             ENDDO    ! ic
-          ENDIF
-       ENDDO    ! ib
-    ENDDO    ! m
+          ENDDO  ! ib
+
+       ENDDO
+
+    ELSE  ! default surfaces
+
+       DO  m = surf_s, surf_e
+
+          k = surf%k(m)
+          norm_fac = 1.0_wp
+
+          IF ( norm )  norm_fac = rho_air_zw(k)  ! normalise vertical fluxes by air density
+
+          DO  ib = 1, nbins_aerosol
+             IF ( aerosol_number(ib)%conc(k,j,i) <= nclim  .OR.  schmidt_num(k+1,ib) < 1.0_wp )    &
+                CYCLE
+
+             SELECT CASE ( depo_surf_par_num )
+
+                CASE ( 1 )
+                   CALL depo_vel_Z01( vc(k+1,ib), surf%us(m), schmidt_num(k+1,ib),                 &
+                                      ra_dry(k,j,i,ib), alpha, gamma, par_a, depo(ib) )
+                CASE ( 2 )
+                   CALL depo_vel_P10( vc(k+1,ib), mag_u(k+1), surf%us(m), kvis(k+1),               &
+                                      schmidt_num(k+1,ib), ra_dry(k,j,i,ib), par_l,                &
+                                      c_brownian_diff, c_interception, c_impaction, beta_im,       &
+                                      c_turb_impaction, depo(ib) )
+             END SELECT
+!
+!--          Calculate changes in surface fluxes due to dry deposition
+             IF ( include_emission )  THEN
+                surf%answs(m,ib) = aerosol_number(ib)%source(j,i) - MAX( 0.0_wp,                   &
+                                   depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i) )
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   surf%amsws(m,icc) = aerosol_mass(icc)%source(j,i) - MAX( 0.0_wp,                &
+                                       depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i) )
+                ENDDO  ! ic
+             ELSE
+                surf%answs(m,ib) = -depo(ib) * norm_fac * aerosol_number(ib)%conc(k,j,i)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   surf%amsws(m,icc) = -depo(ib) *  norm_fac * aerosol_mass(icc)%conc(k,j,i)
+                ENDDO  ! ic
+             ENDIF
+          ENDDO  ! ib
+       ENDDO
+
+    ENDIF
 
  END SUBROUTINE depo_surf
@@ -7200 +7421 @@
  END SUBROUTINE salsa_actions_ij
 
-!
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Call for all grid points
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_non_advective_processes
+
+    USE cpulog,                                                                                    &
+        ONLY:  cpu_log, log_point_s
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i  !<
+    INTEGER(iwp) ::  j  !<
+
+    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
+       IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  THEN
+!
+!--       Calculate aerosol dynamic processes. salsa_driver can be run with a longer time step.
+          CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                CALL salsa_diagnostics( i, j )
+                CALL salsa_driver( i, j, 3 )
+                CALL salsa_diagnostics( i, j )
+             ENDDO
+          ENDDO
+          CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE salsa_non_advective_processes
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Call for grid points i,j
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_non_advective_processes_ij( i, j )
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point_s
+
+    IMPLICIT NONE
+
+    INTEGER(iwp), INTENT(IN) ::  i  !< grid index in x-direction
+    INTEGER(iwp), INTENT(IN) ::  j  !< grid index in y-direction
+
+    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
+       IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  THEN
+!
+!--       Calculate aerosol dynamic processes. salsa_driver can be run with a longer time step.
+          CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
+          CALL salsa_diagnostics( i, j )
+          CALL salsa_driver( i, j, 3 )
+          CALL salsa_diagnostics( i, j )
+          CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE salsa_non_advective_processes_ij
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine for exchange horiz of salsa variables.
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_exchange_horiz_bounds
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point_s
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  ib   !<
+    INTEGER(iwp) ::  ic   !<
+    INTEGER(iwp) ::  icc  !<
+    INTEGER(iwp) ::  ig   !<
+
+    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
+       IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  THEN
+
+          CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
+!
+!--       Exchange ghost points and decycle if needed.
+          DO  ib = 1, nbins_aerosol
+             CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
+             CALL salsa_boundary_conds( aerosol_number(ib)%conc, aerosol_number(ib)%init )
+             DO  ic = 1, ncomponents_mass
+                icc = ( ic - 1 ) * nbins_aerosol + ib
+                CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
+                CALL salsa_boundary_conds( aerosol_mass(icc)%conc, aerosol_mass(icc)%init )
+             ENDDO
+          ENDDO
+          IF ( .NOT. salsa_gases_from_chem )  THEN
+             DO  ig = 1, ngases_salsa
+                CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
+                CALL salsa_boundary_conds( salsa_gas(ig)%conc, salsa_gas(ig)%init )
+             ENDDO
+          ENDIF
+          CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
+!
+!--       Update last_salsa_time
+          last_salsa_time = time_since_reference_point
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE salsa_exchange_horiz_bounds
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -7272 +7602 @@
 ! ------------
 !> Calculate the prognostic equation for aerosol number and mass, and gas
-!> concentrations. Cache-optimized.
+!> concentrations. For vector machines.
 !------------------------------------------------------------------------------!
  SUBROUTINE salsa_prognostic_equations()
@@ -7379 +7709 @@
     REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  flux_l  !<
 
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  rs_p    !<
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  rs      !<
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  trs_m   !<
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  rs_p    !<
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  rs      !<
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  trs_m   !<
 
     icc = ( ic - 1 ) * nbins_aerosol + ib
@@ -7433 +7763 @@
     END SELECT
 !
-!-- Sedimentation for aerosol number and mass
+!-- Sedimentation and prognostic equation for aerosol number and mass
     IF ( lsdepo  .AND.  do_sedimentation )  THEN
-       tend(nzb+1:nzt,j,i) = tend(nzb+1:nzt,j,i) - MAX( 0.0_wp,                                    &
-                             ( rs(nzb+2:nzt+1,j,i) * sedim_vd(nzb+2:nzt+1,j,i,ib) -                &
-                               rs(nzb+1:nzt,j,i) * sedim_vd(nzb+1:nzt,j,i,ib) ) * ddzu(nzb+1:nzt) )&
-                             * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(nzb:nzt-1,j,i), 0 ) )
+!DIR$ IVDEP
+       DO  k = nzb+1, nzt
+          tend(k,j,i) = tend(k,j,i) - MAX( 0.0_wp, ( rs(k+1,j,i) * sedim_vd(k+1,j,i,ib) -          &
+                                                     rs(k,j,i) * sedim_vd(k,j,i,ib) ) * ddzu(k) )  &
+                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k-1,j,i), 0 ) )
+          rs_p(k,j,i) = rs(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * trs_m(k,j,i) )     &
+                                      - tsc(5) * rdf_sc(k) * ( rs(k,j,i) - rs_init(k) ) )          &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+          IF ( rs_p(k,j,i) < 0.0_wp )  rs_p(k,j,i) = 0.1_wp * rs(k,j,i)
+       ENDDO
+    ELSE
+!
+!--    Prognostic equation
+!DIR$ IVDEP
+       DO  k = nzb+1, nzt
+          rs_p(k,j,i) = rs(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * trs_m(k,j,i) )     &
+                                                - tsc(5) * rdf_sc(k) * ( rs(k,j,i) - rs_init(k) ) )&
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+          IF ( rs_p(k,j,i) < 0.0_wp )  rs_p(k,j,i) = 0.1_wp * rs(k,j,i)
+       ENDDO
     ENDIF
-!
-!-- Prognostic equation for a scalar
-    DO  k = nzb+1, nzt
-       rs_p(k,j,i) = rs(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * trs_m(k,j,i) )        &
-                                           - tsc(5) * rdf_sc(k) * ( rs(k,j,i) - rs_init(k) ) )     &
-                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-       IF ( rs_p(k,j,i) < 0.0_wp )  rs_p(k,j,i) = 0.1_wp * rs(k,j,i)
-    ENDDO
 !
 !-- Calculate tendencies for the next Runge-Kutta step
@@ -7468 +7806 @@
 ! ------------
 !> Calculate the tendencies for aerosol number and mass concentrations.
-!> Vector-optimized.
+!> For vector machines.
 !------------------------------------------------------------------------------!
  SUBROUTINE salsa_tendency( id, rs_p, rs, trs_m, ib, ic, rs_init, do_sedimentation )
@@ -9084 +9422 @@
     SELECT CASE ( TRIM( var ) )
 
-       CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV', 'N_bin1', 'N_bin2', 'N_bin3',     &
-              'N_bin4',  'N_bin5',  'N_bin6', 'N_bin7', 'N_bin8', 'N_bin9', 'N_bin10', 'N_bin11',  &
-              'N_bin12', 'Ntot' )
+       CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
           IF (  .NOT.  salsa )  THEN
              message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
-             CALL message( 'check_parameters', 'PA0645', 1, 2, 0, 6, 0 )
+             CALL message( 'check_parameters', 'PA0652', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF (  salsa_gases_from_chem )  THEN
+             message_string = 'gases are imported from the chemistry module and thus output of "'  &
+                               // TRIM( var ) // '" is not allowed'
+             CALL message( 'check_parameters', 'PA0653', 1, 2, 0, 6, 0 )
           ENDIF
           unit = '#/m3'
@@ -9108 +9449 @@
           ENDIF
           unit = 'kg/m3'
+
+       CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4',  'N_bin5',  'N_bin6', 'N_bin7', 'N_bin8',    &
+              'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12', 'Ntot' )
+          IF (  .NOT.  salsa )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
+             CALL message( 'check_parameters', 'PA0645', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = '#/m3'
 
        CASE DEFAULT