Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4268)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4270)
@@ -26,4 +26,15 @@
 ! -----------------
 ! $Id$
+! - Implement offline nesting for salsa
+! - Alphabetic ordering for module interfaces
+! - Remove init_aerosol_type and init_gases_type from salsa_parin and define them
+!   based on the initializing_actions
+! - parameter definition removed from "season" and "season_z01" is added to parin
+! - bugfix in application of index_hh after implementing the new
+!   palm_date_time_mod
+! - Reformat salsa emission data with LOD=2: size distribution given for each
+!   emission category
+! 
+! 4268 2019-10-17 11:29:38Z schwenkel
 ! Moving module specific boundary conditions from time_integration to module
 ! 
@@ -31,5 +42,5 @@
 ! Document previous changes: use global variables nx, ny and nz in salsa_header
 ! 
-! 4255 2019-10-04 11:50:55Z monakurppa
+! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
 ! 
@@ -241,11 +252,8 @@
                                                    !< 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
-                                             !< distribution
+    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 ::  maxspec = 7  !< Max. number of aerosol species
-    INTEGER(iwp), PARAMETER ::  season = 1   !< For dry depostion by Zhang et al.: 1 = summer,
-                                             !< 2 = autumn (no harvest yet), 3 = late autumn
-                                             !< (already frost), 4 = winter, 5 = transitional spring
+
 
     REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
@@ -253,22 +261,16 @@
 !-- Universal constants
     REAL(wp), PARAMETER ::  abo    = 1.380662E-23_wp   !< Boltzmann constant (J/K)
-    REAL(wp), PARAMETER ::  alv    = 2.260E+6_wp       !< latent heat for H2O
-                                                       !< vaporisation (J/kg)
+    REAL(wp), PARAMETER ::  alv    = 2.260E+6_wp       !< latent heat for H2O vaporisation (J/kg)
     REAL(wp), PARAMETER ::  alv_d_rv  = 4896.96865_wp  !< alv / rv
-    REAL(wp), PARAMETER ::  am_airmol = 4.8096E-26_wp  !< Average mass of one air
-                                                       !< molecule (Jacobson,
-                                                       !< 2005, Eq. 2.3)
+    REAL(wp), PARAMETER ::  am_airmol = 4.8096E-26_wp  !< Average mass of an air molecule (Jacobson 2005, Eq.2.3)
     REAL(wp), PARAMETER ::  api6   = 0.5235988_wp      !< pi / 6
     REAL(wp), PARAMETER ::  argas  = 8.314409_wp       !< Gas constant (J/(mol K))
     REAL(wp), PARAMETER ::  argas_d_cpd = 8.281283865E-3_wp  !< argas per cpd
     REAL(wp), PARAMETER ::  avo    = 6.02214E+23_wp    !< Avogadro constant (1/mol)
-    REAL(wp), PARAMETER ::  d_sa   = 5.539376964394570E-10_wp  !< diameter of condensing sulphuric
-                                                               !< acid molecule (m)
+    REAL(wp), PARAMETER ::  d_sa   = 5.539376964394570E-10_wp  !< diameter of condensing H2SO4 molecule (m)
     REAL(wp), PARAMETER ::  for_ppm_to_nconc =  7.243016311E+16_wp !< ppm * avo / R (K/(Pa*m3))
-    REAL(wp), PARAMETER ::  epsoc  = 0.15_wp          !< water uptake of organic
-                                                      !< material
+    REAL(wp), PARAMETER ::  epsoc  = 0.15_wp          !< water uptake of organic material
     REAL(wp), PARAMETER ::  mclim  = 1.0E-23_wp       !< mass concentration min limit (kg/m3)
-    REAL(wp), PARAMETER ::  n3     = 158.79_wp        !< Number of H2SO4 molecules in 3 nm cluster
-                                                      !< if d_sa=5.54e-10m
+    REAL(wp), PARAMETER ::  n3     = 158.79_wp        !< Number of H2SO4 molecules in 3 nm cluster if d_sa=5.54e-10m
     REAL(wp), PARAMETER ::  nclim  = 1.0_wp           !< number concentration min limit (#/m3)
     REAL(wp), PARAMETER ::  surfw0 = 0.073_wp         !< surface tension of water at 293 K (J/m2)
@@ -377,5 +379,4 @@
     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
     INTEGER(iwp) ::  end_subrange_2a = 1    !< last index for bin subrange 2a
@@ -392,12 +393,11 @@
     INTEGER(iwp) ::  init_aerosol_type = 0  !< Initial size distribution type
                                             !< 0 = uniform (read from PARIN)
-                                            !< 1 = read vertical profile of the mode number
-                                            !<     concentration from an input file
+                                            !< 1 = read vertical profiles from an input file
     INTEGER(iwp) ::  init_gases_type = 0    !< Initial gas concentration type
                                             !< 0 = uniform (read from PARIN)
-                                            !< 1 = read vertical profile from an input file
+                                            !< 1 = read vertical profiles from an input file
     INTEGER(iwp) ::  lod_gas_emissions = 0  !< level of detail of the gaseous emission data
-    INTEGER(iwp) ::  main_street_id = 0     !< lower bound of main street IDs (OpenStreetMaps) for parameterized mode
-    INTEGER(iwp) ::  max_street_id = 0      !< upper bound of main street IDs (OpenStreetMaps) for parameterized mode
+    INTEGER(iwp) ::  main_street_id = 0     !< lower bound of main street IDs for parameterized emission mode
+    INTEGER(iwp) ::  max_street_id = 0      !< upper bound of main street IDs for parameterized emission mode
     INTEGER(iwp) ::  nbins_aerosol = 1      !< total number of size bins
     INTEGER(iwp) ::  ncc   = 1              !< number of chemical components used
@@ -422,5 +422,8 @@
                                             !<     + heteromolecular nucleation with H2SO4*ORG
     INTEGER(iwp) ::  salsa_pr_count = 0     !< counter for salsa variable profiles
-    INTEGER(iwp) ::  side_street_id = 0     !< lower bound of side street IDs (OpenStreetMaps) for parameterized mode
+    INTEGER(iwp) ::  season_z01 = 1         !< For dry deposition by Zhang et al.: 1 = summer,
+                                            !< 2 = autumn (no harvest yet), 3 = late autumn
+                                            !< (already frost), 4 = winter, 5 = transitional spring
+    INTEGER(iwp) ::  side_street_id = 0     !< lower bound of side street IDs for parameterized emission mode
     INTEGER(iwp) ::  start_subrange_1a = 1  !< start index for bin subranges: subrange 1a
     INTEGER(iwp) ::  start_subrange_2a = 1  !<                                subrange 2a
@@ -448,10 +451,9 @@
     LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
     LOGICAL ::  nest_salsa              = .FALSE.  !< Apply nesting for salsa
+    LOGICAL ::  nesting_offline_salsa   = .FALSE.  !< Apply offline 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 ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA
+    LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Include van der Waals and viscous forces in coagulation
     LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
 !
@@ -491,4 +493,5 @@
     REAL(wp) ::  ocsv_init = nclim                   !< Init value for semi-volatile organic gases
     REAL(wp) ::  rhlim = 1.20_wp                     !< RH limit in %/100. Prevents unrealistical RH
+    REAL(wp) ::  time_utc_init                       !< time in seconds-of-day of origin_date_time
     REAL(wp) ::  skip_time_do_salsa = 0.0_wp         !< Starting time of SALSA (s)
 !
@@ -503,9 +506,9 @@
 !
 !-- Initial mass fractions / chemical composition of the size distribution
-    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_a = & !< mass fractions between
-             (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for A bins
-    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_b = & !< mass fractions between
-             (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for B bins
-    REAL(wp), DIMENSION(nreg+1) ::  reglim = & !< Min&max diameters of size subranges
+    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_a = &  !< mass fractions between
+             (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)     !< aerosol species for A bins
+    REAL(wp), DIMENSION(maxspec) ::  mass_fracs_b = &  !< mass fractions between
+             (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)     !< aerosol species for B bins
+    REAL(wp), DIMENSION(nreg+1) ::  reglim = &         !< Min&max diameters of size subranges
                                  (/ 3.0E-9_wp, 5.0E-8_wp, 1.0E-5_wp/)
 !
@@ -628,12 +631,58 @@
        REAL(wp) ::  fill  !< fill value
 
-       REAL(wp), DIMENSION(:), ALLOCATABLE :: preproc_mass_fracs  !< mass fractions
-
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: def_mass_fracs  !< mass fractions per emis. category
-
-       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: def_data      !< surface emission values in PM
-       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: preproc_data  !< surface emission values per bin
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mass_fracs  !< mass fractions per emis. category
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: num_fracs   !< number fractions per emis. category
+
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: def_data      !< surface emission in PM
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: preproc_data  !< surface emission per category
 
     END TYPE salsa_emission_value_type
+!
+!-- Offline nesting data type
+    TYPE salsa_nest_offl_type
+
+       CHARACTER(LEN=16) ::  char_l = 'ls_forcing_left_'  !< leading substring at left boundary
+       CHARACTER(LEN=17) ::  char_n = 'ls_forcing_north_' !< leading substring at north boundary
+       CHARACTER(LEN=17) ::  char_r = 'ls_forcing_right_' !< leading substring at right boundary
+       CHARACTER(LEN=17) ::  char_s = 'ls_forcing_south_' !< leading substring at south boundary
+       CHARACTER(LEN=15) ::  char_t = 'ls_forcing_top_'   !< leading substring at top boundary
+
+       CHARACTER(LEN=5), DIMENSION(1:ngases_salsa) ::  gas_name = (/'h2so4','hno3 ','nh3  ','ocnv ','ocsv '/)
+
+       CHARACTER(LEN=25),  DIMENSION(:), ALLOCATABLE ::  cc_name    !< chemical component name
+       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< list of variable names
+
+       INTEGER(iwp) ::  id_dynamic  !< NetCDF id of dynamic input file
+       INTEGER(iwp) ::  ncc         !< number of aerosol chemical components
+       INTEGER(iwp) ::  nt          !< number of time levels in dynamic input file
+       INTEGER(iwp) ::  nzu         !< number of vertical levels on scalar grid in dynamic input file
+       INTEGER(iwp) ::  tind        !< time index for reference time in mesoscale-offline nesting
+       INTEGER(iwp) ::  tind_p      !< time index for following time in mesoscale-offline nesting
+
+       INTEGER(iwp), DIMENSION(maxspec) ::  cc_in2mod = 0  !< to transfer chemical composition from input to model
+
+       LOGICAL ::  init  = .FALSE. !< flag indicating the initialisation of offline nesting
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dmid      !< vertical profile of aerosol bin diameters
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time      !< time in dynamic input file
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu_atmos  !< zu in dynamic input file
+
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  gconc_left   !< gas conc. at left boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  gconc_north  !< gas conc. at north boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  gconc_right  !< gas conc. at right boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  gconc_south  !< gas conc. at south boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  gconc_top    !< gas conc.at top boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mconc_left   !< aerosol mass conc. at left boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mconc_north  !< aerosol mass conc. at north boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mconc_right  !< aerosol mass conc. at right boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mconc_south  !< aerosol mass conc. at south boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mconc_top    !< aerosol mass conc. at top boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  nconc_left   !< aerosol number conc. at left boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  nconc_north  !< aerosol number conc. at north boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  nconc_right  !< aerosol number conc. at right boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  nconc_south  !< aerosol number conc. at south boundary
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  nconc_top    !< aerosol number conc. at top boundary
+
+    END TYPE salsa_nest_offl_type
 !
 !-- Prognostic variable: Aerosol size bin information (number (#/m3) and mass (kg/m3) concentration)
@@ -728,4 +777,7 @@
                                     !< component name: 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
 !
+!-- Offline nesting:
+    TYPE(salsa_nest_offl_type) ::  salsa_nest_offl  !< data structure for offline nesting
+!
 !-- Data output arrays:
 !
@@ -808,8 +860,4 @@
     END INTERFACE salsa_define_netcdf_grid
 
-    INTERFACE salsa_driver
-       MODULE PROCEDURE salsa_driver
-    END INTERFACE salsa_driver
-
     INTERFACE salsa_emission_update
        MODULE PROCEDURE salsa_emission_update
@@ -831,4 +879,16 @@
        MODULE PROCEDURE salsa_init_arrays
     END INTERFACE salsa_init_arrays
+
+    INTERFACE salsa_nesting_offl_bc
+       MODULE PROCEDURE salsa_nesting_offl_bc
+    END INTERFACE salsa_nesting_offl_bc
+
+    INTERFACE salsa_nesting_offl_init
+       MODULE PROCEDURE salsa_nesting_offl_init
+    END INTERFACE salsa_nesting_offl_init
+
+    INTERFACE salsa_nesting_offl_input
+       MODULE PROCEDURE salsa_nesting_offl_input
+    END INTERFACE salsa_nesting_offl_input
 
     INTERFACE salsa_non_advective_processes
@@ -873,20 +933,56 @@
 !
 !-- Public functions:
-    PUBLIC salsa_boundary_conds, salsa_check_data_output, salsa_check_parameters, salsa_boundary_conditions, &
-           salsa_3d_data_averaging, salsa_data_output_2d, salsa_data_output_3d,                    &
-           salsa_data_output_mask, salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver,      &
-           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_non_advective_processes, salsa_exchange_horiz_bounds,              &
-           salsa_check_data_output_pr, salsa_statistics
+    PUBLIC salsa_3d_data_averaging,       &
+           salsa_actions,                 &
+           salsa_boundary_conds,          &
+           salsa_boundary_conditions,     &
+           salsa_check_data_output,       &
+           salsa_check_data_output_pr,    &
+           salsa_check_parameters,        &
+           salsa_data_output_2d,          &
+           salsa_data_output_3d,          &
+           salsa_data_output_mask,        &
+           salsa_define_netcdf_grid,      &
+           salsa_diagnostics,             &
+           salsa_emission_update,         &
+           salsa_exchange_horiz_bounds,   &
+           salsa_header,                  &
+           salsa_init,                    &
+           salsa_init_arrays,             &
+           salsa_nesting_offl_bc,         &
+           salsa_nesting_offl_init,       &
+           salsa_nesting_offl_input,      &
+           salsa_non_advective_processes, &
+           salsa_parin,                   &
+           salsa_prognostic_equations,    &
+           salsa_rrd_local,               &
+           salsa_statistics,              &
+           salsa_swap_timelevel,          &
+           salsa_wrd_local
+
 !
 !-- Public parameters, constants and initial values
-    PUBLIC bc_am_t_val, bc_an_t_val, bc_gt_t_val, dots_salsa, dt_salsa,                            &
-           ibc_salsa_b, last_salsa_time, lsdepo, nest_salsa, salsa, salsa_gases_from_chem,         &
+    PUBLIC bc_am_t_val,           &
+           bc_an_t_val,           &
+           bc_gt_t_val,           &
+           ibc_salsa_b,           &
+           init_aerosol_type,     &
+           init_gases_type,       &
+           nest_salsa,            &
+           nesting_offline_salsa, &
+           salsa_gases_from_chem, &
            skip_time_do_salsa
 !
-!-- Public prognostic variables
-    PUBLIC aerosol_mass, aerosol_number, gconc_2, mconc_2, nbins_aerosol, ncc, ncomponents_mass,   &
-           nclim, nconc_2, ngases_salsa, prtcl, ra_dry, salsa_gas, sedim_vd
+!-- Public variables
+    PUBLIC aerosol_mass,     &
+           aerosol_number,   &
+           gconc_2,          &
+           mconc_2,          &
+           nbins_aerosol,    &
+           ncomponents_mass, &
+           nconc_2,          &
+           ngases_salsa,     &
+           salsa_gas,        &
+           salsa_nest_offl
 
 
@@ -930,6 +1026,4 @@
                                      h2so4_init,                               &
                                      hno3_init,                                &
-                                     init_gases_type,                          &
-                                     init_aerosol_type,                        &
                                      listspec,                                 &
                                      main_street_id,                           &
@@ -940,4 +1034,5 @@
                                      nbin,                                     &
                                      nest_salsa,                               &
+                                     nesting_offline_salsa,                    &
                                      nf2a,                                     &
                                      nh3_init,                                 &
@@ -958,4 +1053,5 @@
                                      salsa,                                    &
                                      salsa_emission_mode,                      &
+                                     season_z01,                               &
                                      sigmag,                                   &
                                      side_street_id,                           &
@@ -1002,15 +1098,16 @@
 
     USE control_parameters,                                                                        &
-        ONLY:  humidity
+        ONLY:  humidity, initializing_actions
 
     IMPLICIT NONE
 
 !
-!-- Checks go here (cf. check_parameters.f90).
+!-- Check that humidity is switched on
     IF ( salsa  .AND.  .NOT.  humidity )  THEN
        WRITE( message_string, * ) 'salsa = ', salsa, ' is not allowed with humidity = ', humidity
        CALL message( 'salsa_check_parameters', 'PA0594', 1, 2, 0, 6, 0 )
     ENDIF
-
+!
+!-- Set bottom boundary condition flag
     IF ( bc_salsa_b == 'dirichlet' )  THEN
        ibc_salsa_b = 0
@@ -1021,5 +1118,6 @@
        CALL message( 'salsa_check_parameters', 'PA0595', 1, 2, 0, 6, 0 )
     ENDIF
-
+!
+!-- Set top boundary conditions flag
     IF ( bc_salsa_t == 'dirichlet' )  THEN
        ibc_salsa_t = 0
@@ -1032,16 +1130,32 @@
        CALL message( 'salsa_check_parameters', 'PA0596', 1, 2, 0, 6, 0 )
     ENDIF
-
+!
+!-- If nest_salsa = .F., set top boundary to dirichlet
+    IF ( .NOT. nest_salsa  .AND.  ibc_salsa_t == 2  )  THEN
+       ibc_salsa_t = 0
+       bc_salsa_t = 'dirichlet'
+    ENDIF
+!
+!-- Check J3 parametrisation
     IF ( nj3 < 1  .OR.  nj3 > 3 )  THEN
        message_string = 'unknown nj3 (must be 1-3)'
        CALL message( 'salsa_check_parameters', 'PA0597', 1, 2, 0, 6, 0 )
     ENDIF
-
+!
+!-- Check bottom boundary condition in case of surface emissions
     IF ( salsa_emission_mode /= 'no_emission'  .AND.  ibc_salsa_b  == 0 ) THEN
        message_string = 'salsa_emission_mode /= "no_emission" requires bc_salsa_b = "Neumann"'
        CALL message( 'salsa_check_parameters','PA0598', 1, 2, 0, 6, 0 )
     ENDIF
-
+!
+!-- Check whether emissions are applied
     IF ( salsa_emission_mode /= 'no_emission' )  include_emission = .TRUE.
+!
+!-- Set the initialisation type: background concentration are read from PIDS_DYNAMIC if
+!-- initializing_actions = 'inifor set_constant_profiles'
+    IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
+       init_aerosol_type = 1
+       init_gases_type = 1
+    ENDIF
 
  END SUBROUTINE salsa_check_parameters
@@ -3751,5 +3865,5 @@
        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
+       par_a   = A_z01(depo_pcm_type_num, season_z01) * 1.0E-3_wp
 !
 !--    Deposition efficiencies from Table 1. Constants from Table 2.
@@ -3992,5 +4106,5 @@
     alpha   = alpha_z01(luc_urban)
     gamma   = gamma_z01(luc_urban)
-    par_a   = A_z01(luc_urban, season) * 1.0E-3_wp
+    par_a   = A_z01(luc_urban, season_z01) * 1.0E-3_wp
 
     par_l            = l_p10(luc_urban) * 0.01_wp
@@ -4014,5 +4128,5 @@
              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
+             par_a = A_z01( match_array%match_lupg(m), season_z01 ) * 1.0E-3_wp
 
              beta_im          = beta_im_p10( match_array%match_lupg(m) )
@@ -4045,5 +4159,5 @@
              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
+             par_a = A_z01( match_array%match_luvw(m), season_z01 ) * 1.0E-3_wp
 
              beta_im          = beta_im_p10( match_array%match_luvw(m) )
@@ -4076,5 +4190,5 @@
              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
+             par_a = A_z01( match_array%match_luww(m), season_z01 ) * 1.0E-3_wp
 
              beta_im          = beta_im_p10( match_array%match_luww(m) )
@@ -8215,4 +8329,7 @@
  SUBROUTINE salsa_boundary_conds_decycle ( sq, sq_init )
 
+    USE control_parameters,                                                                        &
+        ONLY:  nesting_offline
+
     IMPLICIT NONE
 
@@ -8232,4 +8349,7 @@
 
     flag = 0.0_wp
+!
+!-- Skip input if forcing from larger-scale models is applied.
+    IF ( nesting_offline  .AND.  nesting_offline_salsa )  RETURN
 !
 !-- Left and right boundaries
@@ -8403,4 +8523,7 @@
  SUBROUTINE salsa_emission_update
 
+    USE palm_date_time_mod,                                                                        &
+        ONLY:  get_date_time
+
     IMPLICIT NONE
 
@@ -8408,10 +8531,15 @@
 
        IF ( time_since_reference_point >= skip_time_do_salsa  )  THEN
-
-          IF ( next_aero_emission_update <= time_since_reference_point )  THEN
+!
+!--       Get time_utc_init from origin_date_time
+          CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
+
+          IF ( next_aero_emission_update <=                                                        &
+               MAX( time_since_reference_point, 0.0_wp ) + time_utc_init )  THEN
              CALL salsa_emission_setup( .FALSE. )
           ENDIF
 
-          IF ( next_gas_emission_update <= time_since_reference_point )  THEN
+          IF ( next_gas_emission_update <=                                                         &
+               MAX( time_since_reference_point, 0.0_wp ) + time_utc_init )  THEN
              IF ( salsa_emission_mode == 'read_from_file'  .AND.  .NOT. salsa_gases_from_chem )    &
              THEN
@@ -8474,5 +8602,11 @@
     LOGICAL, INTENT(in) ::  init  !< if .TRUE. --> initialisation call
 
-    REAL(wp) ::  second_of_day    !< second of the day
+    REAL(wp) ::  second_of_day  !< second of the day
+
+    REAL(wp), DIMENSION(24) ::  par_emis_time_factor =  & !< time factors for the parameterized mode
+                                                      (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, &
+                                                         0.056, 0.053, 0.051, 0.051, 0.052, 0.055, &
+                                                         0.059, 0.061, 0.064, 0.067, 0.069, 0.069, &
+                                                         0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  nsect_emission  !< sectional number emission
@@ -8501,4 +8635,12 @@
                 ENDDO
              ELSE
+!
+!--             Get a time factor for the specific hour
+                IF ( .NOT.  ALLOCATED( aero_emission_att%time_factor ) )                           &
+                   ALLOCATE( aero_emission_att%time_factor(1) )
+                CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), hour=hour_of_day )
+                index_hh = hour_of_day
+                aero_emission_att%time_factor(1) = par_emis_time_factor(index_hh+1)
+
                 IF ( street_type_f%from_file )  THEN
                    DO  i = nxl, nxr
@@ -8506,8 +8648,10 @@
                          IF ( street_type_f%var(j,i) >= main_street_id  .AND.                      &
                               street_type_f%var(j,i) < max_street_id )  THEN
-                            source_array(j,i,:) = nsect_emission(:) * emiss_factor_main
+                            source_array(j,i,:) = nsect_emission(:) * emiss_factor_main *          &
+                                                  aero_emission_att%time_factor(1)
                          ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.                  &
                                   street_type_f%var(j,i) < main_street_id )  THEN
-                            source_array(j,i,:) = nsect_emission(:) * emiss_factor_side
+                            source_array(j,i,:) = nsect_emission(:) * emiss_factor_side *          &
+                                                  aero_emission_att%time_factor(1)
                          ENDIF
                       ENDDO
@@ -8587,6 +8731,5 @@
 !
 !--          Read the index and name of chemical components
-             CALL get_dimension_length( id_salsa, aero_emission_att%ncc,         &
-                                                          'composition_index' )
+             CALL get_dimension_length( id_salsa, aero_emission_att%ncc, 'composition_index' )
              ALLOCATE( aero_emission_att%cc_index(1:aero_emission_att%ncc) )
              CALL get_variable( id_salsa, 'composition_index', aero_emission_att%cc_index )
@@ -8629,4 +8772,42 @@
              ENDIF
 !
+!--          Get number of emission categories
+             CALL get_dimension_length( id_salsa, aero_emission_att%ncat, 'ncat' )
+!
+!--          Get the chemical composition (i.e. mass fraction of different species) in aerosols
+             IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) )  THEN
+                ALLOCATE( aero_emission%mass_fracs(1:aero_emission_att%ncat,                       &
+                                                   1:aero_emission_att%ncc) )
+                CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%mass_fracs,      &
+                                   0, aero_emission_att%ncc-1, 0, aero_emission_att%ncat-1 )
+             ELSE
+                message_string = 'Missing emission_mass_fracs in ' //  TRIM( input_file_salsa )
+                CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 )
+             ENDIF
+!
+!--          If the chemical component is not activated, set its mass fraction to 0 to avoid
+!--          inbalance between number and mass flux
+             cc_i2m = aero_emission_att%cc_in2mod
+             IF ( index_so4 < 0  .AND.  cc_i2m(1) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(1)) = 0.0_wp
+             IF ( index_oc  < 0  .AND.  cc_i2m(2) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(2)) = 0.0_wp
+             IF ( index_bc  < 0  .AND.  cc_i2m(3) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(3)) = 0.0_wp
+             IF ( index_du  < 0  .AND.  cc_i2m(4) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(4)) = 0.0_wp
+             IF ( index_ss  < 0  .AND.  cc_i2m(5) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(5)) = 0.0_wp
+             IF ( index_no  < 0  .AND.  cc_i2m(6) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(6)) = 0.0_wp
+             IF ( index_nh  < 0  .AND.  cc_i2m(7) > 0 )                                            &
+                aero_emission%mass_fracs(:,cc_i2m(7)) = 0.0_wp
+!
+!--          Then normalise the mass fraction so that SUM = 1
+             DO  in = 1, aero_emission_att%ncat
+                aero_emission%mass_fracs(in,:) = aero_emission%mass_fracs(in,:) /                  &
+                                                 SUM( aero_emission%mass_fracs(in,:) )
+             ENDDO
+!
 !--          Inquire the fill value
              CALL get_attribute( id_salsa, '_FillValue', aero_emission%fill, .FALSE.,              &
@@ -8659,6 +8840,5 @@
                 ENDIF
 !
-!--             Get number of emission categories and allocate emission arrays
-                CALL get_dimension_length( id_salsa, aero_emission_att%ncat, 'ncat' )
+!--             Allocate emission arrays
                 ALLOCATE( aero_emission_att%cat_index(1:aero_emission_att%ncat),                   &
                           aero_emission_att%rho(1:aero_emission_att%ncat),                         &
@@ -8697,6 +8877,5 @@
 !--             For each hour of year:
                 IF ( check_existence( aero_emission_att%var_names, 'nhoursyear' ) )  THEN
-                   CALL get_dimension_length( id_salsa,                                            &
-                                              aero_emission_att%nhoursyear, 'nhoursyear' )
+                   CALL get_dimension_length( id_salsa, aero_emission_att%nhoursyear, 'nhoursyear' )
                    ALLOCATE( aero_emission_att%etf(1:aero_emission_att%ncat,                       &
                                                    1:aero_emission_att%nhoursyear) )
@@ -8706,6 +8885,5 @@
 !--             Based on the month, day and hour:
                 ELSEIF ( check_existence( aero_emission_att%var_names, 'nmonthdayhour' ) )  THEN
-                   CALL get_dimension_length( id_salsa,                                            &
-                                              aero_emission_att%nmonthdayhour,                     &
+                   CALL get_dimension_length( id_salsa, aero_emission_att%nmonthdayhour,           &
                                               'nmonthdayhour' )
                    ALLOCATE( aero_emission_att%etf(1:aero_emission_att%ncat,                       &
@@ -8720,42 +8898,6 @@
 !
 !--             Next emission update
-                CALL get_date_time( 0.0_wp, second_of_day=second_of_day )
-                next_aero_emission_update = MOD( second_of_day, seconds_per_hour ) - seconds_per_hour
-!
-!--             Get chemical composition (i.e. mass fraction of different species) in aerosols
-                IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) )  THEN
-                   ALLOCATE( aero_emission%def_mass_fracs(1:aero_emission_att%ncat,                &
-                                                          1:aero_emission_att%ncc) )
-                   aero_emission%def_mass_fracs = 0.0_wp
-                   CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%def_mass_fracs,&
-                                      0, aero_emission_att%ncc-1, 0, aero_emission_att%ncat-1 )
-                ELSE
-                   message_string = 'Missing emission_mass_fracs in ' //  TRIM( input_file_salsa )
-                   CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 )
-                ENDIF
-!
-!--             If the chemical component is not activated, set its mass fraction to 0 to avoid
-!--             inbalance between number and mass flux
-                cc_i2m = aero_emission_att%cc_in2mod
-                IF ( index_so4 < 0  .AND.  cc_i2m(1) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(1)) = 0.0_wp
-                IF ( index_oc  < 0  .AND.  cc_i2m(2) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(2)) = 0.0_wp
-                IF ( index_bc  < 0  .AND.  cc_i2m(3) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(3)) = 0.0_wp
-                IF ( index_du  < 0  .AND.  cc_i2m(4) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(4)) = 0.0_wp
-                IF ( index_ss  < 0  .AND.  cc_i2m(5) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(5)) = 0.0_wp
-                IF ( index_no  < 0  .AND.  cc_i2m(6) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(6)) = 0.0_wp
-                IF ( index_nh  < 0  .AND.  cc_i2m(7) > 0 )                                         &
-                                                  aero_emission%def_mass_fracs(:,cc_i2m(7)) = 0.0_wp
-!
-!--             Then normalise the mass fraction so that SUM = 1
-                DO  in = 1, aero_emission_att%ncat
-                   aero_emission%def_mass_fracs(in,:) = aero_emission%def_mass_fracs(in,:) /       &
-                                                       SUM( aero_emission%def_mass_fracs(in,:) )
-                ENDDO
+                CALL get_date_time( time_since_reference_point, second_of_day=second_of_day )
+                next_aero_emission_update = MOD( second_of_day, seconds_per_hour ) !- seconds_per_hour
 !
 !--             Calculate average mass density (kg/m3)
@@ -8763,17 +8905,17 @@
 
                 IF ( cc_i2m(1) /= 0 )  aero_emission_att%rho = aero_emission_att%rho +  arhoh2so4 *&
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(1))
+                                                               aero_emission%mass_fracs(:,cc_i2m(1))
                 IF ( cc_i2m(2) /= 0 )  aero_emission_att%rho = aero_emission_att%rho + arhooc *    &
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(2))
+                                                               aero_emission%mass_fracs(:,cc_i2m(2))
                 IF ( cc_i2m(3) /= 0 )  aero_emission_att%rho = aero_emission_att%rho + arhobc *    &
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(3))
+                                                               aero_emission%mass_fracs(:,cc_i2m(3))
                 IF ( cc_i2m(4) /= 0 )  aero_emission_att%rho = aero_emission_att%rho + arhodu *    &
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(4))
+                                                               aero_emission%mass_fracs(:,cc_i2m(4))
                 IF ( cc_i2m(5) /= 0 )  aero_emission_att%rho = aero_emission_att%rho + arhoss *    &
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(5))
+                                                               aero_emission%mass_fracs(:,cc_i2m(5))
                 IF ( cc_i2m(6) /= 0 )  aero_emission_att%rho = aero_emission_att%rho + arhohno3 *  &
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(6))
+                                                               aero_emission%mass_fracs(:,cc_i2m(6))
                 IF ( cc_i2m(7) /= 0 )  aero_emission_att%rho = aero_emission_att%rho + arhonh3 *   &
-                                                           aero_emission%def_mass_fracs(:,cc_i2m(7))
+                                                               aero_emission%mass_fracs(:,cc_i2m(7))
 !
 !--             Allocate and read surface emission data (in total PM)
@@ -8809,5 +8951,5 @@
                 ALLOCATE( aero_emission_att%dmid(1:nbins_aerosol),                                 &
                           aero_emission_att%time(1:aero_emission_att%nt),                          &
-                          aero_emission%preproc_mass_fracs(1:aero_emission_att%ncc) )
+                          aero_emission%num_fracs(1:aero_emission_att%ncat,1:nbins_aerosol) )
 !
 !--             Read mean diameters
@@ -8831,38 +8973,18 @@
                 ENDIF
 !
-!--             Read emission mass fractions
-                IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) )  THEN
-                   CALL get_variable( id_salsa, 'emission_mass_fracs',                             &
-                                      aero_emission%preproc_mass_fracs )
+!--             Read emission number fractions per category
+                IF ( check_existence( aero_emission_att%var_names, 'emission_number_fracs' ) )  THEN
+                   CALL get_variable( id_salsa, 'emission_number_fracs', aero_emission%num_fracs,  &
+                                      0, nbins_aerosol-1, 0, aero_emission_att%ncat-1 )
                 ELSE
-                   message_string = 'Missing emission_mass_fracs in ' //  TRIM( input_file_salsa )
+                   message_string = 'Missing emission_number_fracs in ' //  TRIM( input_file_salsa )
                    CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 )
                 ENDIF
-!
-!--             If the chemical component is not activated, set its mass fraction to 0
-                cc_i2m = aero_emission_att%cc_in2mod
-                IF ( index_so4 < 0  .AND.  cc_i2m(1) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(1)) = 0.0_wp
-                IF ( index_oc  < 0  .AND.  cc_i2m(2) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(2)) = 0.0_wp
-                IF ( index_bc  < 0  .AND.  cc_i2m(3) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(3)) = 0.0_wp
-                IF ( index_du  < 0  .AND.  cc_i2m(4) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(4)) = 0.0_wp
-                IF ( index_ss  < 0  .AND.  cc_i2m(5) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(5)) = 0.0_wp
-                IF ( index_no  < 0  .AND.  cc_i2m(6) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(6)) = 0.0_wp
-                IF ( index_nh  < 0  .AND.  cc_i2m(7) /= 0 )                                        &
-                   aero_emission%preproc_mass_fracs(cc_i2m(7)) = 0.0_wp
-!
-!--             Then normalise the mass fraction so that SUM = 1
-                aero_emission%preproc_mass_fracs = aero_emission%preproc_mass_fracs /              &
-                                                   SUM( aero_emission%preproc_mass_fracs )
 
              ELSE
                 message_string = 'Unknown lod for aerosol_emission_values.'
                 CALL message( 'salsa_emission','PA0637', 1, 2, 0, 6, 0 )
-             ENDIF
+
+             ENDIF  ! lod
 
           ENDIF  ! init
@@ -8877,8 +8999,8 @@
 !
 !--             Get the index of the current hour
-                CALL get_date_time( time_since_reference_point, &
+                CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ),                     &
                                     day_of_year=day_of_year, hour=hour_of_day )
                 index_hh = ( day_of_year - 1_iwp ) * hours_per_day + hour_of_day
-                aero_emission_att%time_factor = aero_emission_att%etf(:,index_hh)
+                aero_emission_att%time_factor = aero_emission_att%etf(:,index_hh+1)
 
              ELSEIF ( aero_emission_att%nhoursyear < aero_emission_att%nmonthdayhour )  THEN
@@ -8886,9 +9008,6 @@
 !--             Get the index of current hour (index_hh) (TODO: Now "workday" is always assumed.
 !--             Needs to be calculated.)
-                CALL get_date_time( time_since_reference_point, &
-                                    month=month_of_year,        &
-                                    day=day_of_month,           &
-                                    hour=hour_of_day,           &
-                                    day_of_week=day_of_week     )
+                CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ), month=month_of_year,&
+                                    day=day_of_month, hour=hour_of_day, day_of_week=day_of_week )
                 index_mm = month_of_year
                 index_dd = months_per_year + day_of_week
@@ -8907,5 +9026,5 @@
                 aero_emission_att%time_factor = aero_emission_att%etf(:,index_mm) *                &
                                                 aero_emission_att%etf(:,index_dd) *                &
-                                                aero_emission_att%etf(:,index_hh)
+                                                aero_emission_att%etf(:,index_hh+1)
              ENDIF
 
@@ -8939,10 +9058,10 @@
                 IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
                    CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod,                      &
-                                  aero_emission%def_mass_fracs(in,:), source_array )
+                                  aero_emission%mass_fracs(in,:), source_array )
                 ELSE
                    CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod,                         &
-                                  aero_emission%def_mass_fracs(in,:), source_array )
+                                  aero_emission%mass_fracs(in,:), source_array )
                    CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod,                         &
-                                  aero_emission%def_mass_fracs(in,:), source_array )
+                                  aero_emission%mass_fracs(in,:), source_array )
                 ENDIF
              ENDDO
@@ -8957,32 +9076,46 @@
           ELSEIF ( aero_emission_att%lod == 2 )  THEN
 !
-!--          Obtain time index for current input starting at 0.
-!--          @todo: At the moment emission data and simulated time correspond to each other.
-             aero_emission_att%tind = MINLOC( ABS( aero_emission_att%time -                        &
-                                                   time_since_reference_point ), DIM = 1 ) - 1
+!--          Get time_utc_init from origin_date_time
+             CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
+!
+!--          Obtain time index for current point in time. Note, the time coordinate in the input
+!--          file is relative to time_utc_init.
+             aero_emission_att%tind = MINLOC( ABS( aero_emission_att%time - (                      &
+                                                   time_utc_init + MAX( time_since_reference_point,&
+                                                                        0.0_wp) ) ), DIM = 1 ) - 1
 !
 !--          Allocate the data input array always before reading in the data and deallocate after
-             ALLOCATE( aero_emission%preproc_data(nys:nyn,nxl:nxr,1:nbins_aerosol) )
+             ALLOCATE( aero_emission%preproc_data(nys:nyn,nxl:nxr,1:aero_emission_att%ncat),       &
+                       source_array(nys:nyn,nxl:nxr,1:nbins_aerosol) )
 !
 !--          Read in the next time step
-             CALL get_variable( id_salsa, 'aerosol_emission_values', aero_emission%preproc_data,&
-                                aero_emission_att%tind, 0, nbins_aerosol-1, nxl, nxr, nys, nyn )
-!
-!--          Set surface fluxes of aerosol number and mass on horizontal surfaces. Set fluxes only
-!--          for either default, land and urban surface.
-             IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
-                CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod,                         &
-                               aero_emission%preproc_mass_fracs, aero_emission%preproc_data )
-             ELSE
-                CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod,                            &
-                               aero_emission%preproc_mass_fracs, aero_emission%preproc_data )
-                CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod,                            &
-                               aero_emission%preproc_mass_fracs, aero_emission%preproc_data )
-             ENDIF
+             CALL get_variable( id_salsa, 'aerosol_emission_values', aero_emission%preproc_data,   &
+                                aero_emission_att%tind, 0, aero_emission_att%ncat-1,               &
+                                nxl, nxr, nys, nyn )
+!
+!--          Calculate the sources per category and set surface fluxes
+             source_array = 0.0_wp
+             DO  in = 1, aero_emission_att%ncat
+                DO  ib = 1, nbins_aerosol
+                   source_array(:,:,ib) = aero_emission%preproc_data(:,:,in) *                     &
+                                          aero_emission%num_fracs(in,ib)
+                ENDDO
+!
+!--             Set fluxes only for either default, land and urban surface.
+                IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
+                   CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod,                      &
+                                  aero_emission%mass_fracs(in,:), source_array )
+                ELSE
+                   CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod,                         &
+                                  aero_emission%mass_fracs(in,:), source_array )
+                   CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod,                         &
+                                  aero_emission%mass_fracs(in,:), source_array )
+                ENDIF
+             ENDDO
 !
 !--          Determine the next emission update
              next_aero_emission_update = aero_emission_att%time(aero_emission_att%tind+2)
 
-             DEALLOCATE( aero_emission%preproc_data )
+             DEALLOCATE( aero_emission%preproc_data, source_array )
 
           ENDIF
@@ -9347,5 +9480,5 @@
 !--       Next emission update
           CALL get_date_time( time_since_reference_point, second_of_day=second_of_day )
-          next_gas_emission_update = MOD( second_of_day, seconds_per_hour ) - seconds_per_hour
+          next_gas_emission_update = MOD( second_of_day, seconds_per_hour ) !- seconds_per_hour
 !
 !--       Allocate and read surface emission data (in total PM) (NOTE that "preprocessed" input data
@@ -9395,5 +9528,7 @@
                               day_of_year=day_of_year, hour=hour_of_day )
           index_hh = ( day_of_year - 1_iwp ) * hours_per_day + hour_of_day
-          time_factor = chem_emission_att%hourly_emis_time_factor(:,index_hh)
+          IF ( .NOT. ALLOCATED( time_factor ) )  ALLOCATE( time_factor(1:chem_emission_att%ncat) )
+          time_factor = 0.0_wp
+          time_factor = chem_emission_att%hourly_emis_time_factor(:,index_hh+1)
 
        ELSEIF ( chem_emission_att%nhoursyear < chem_emission_att%nmonthdayhour )  THEN
@@ -9408,5 +9543,5 @@
           index_mm = month_of_year
           index_dd = months_per_year + day_of_week
-          SELECT CASE(TRIM(daytype))
+          SELECT CASE( TRIM( daytype ) )
 
              CASE ("workday")
@@ -9422,5 +9557,5 @@
           time_factor = chem_emission_att%mdh_emis_time_factor(:,index_mm) *                       &
                         chem_emission_att%mdh_emis_time_factor(:,index_dd) *                       &
-                        chem_emission_att%mdh_emis_time_factor(:,index_hh)
+                        chem_emission_att%mdh_emis_time_factor(:,index_hh+1)
        ENDIF
 !
@@ -9435,11 +9570,11 @@
 !--       Set surface fluxes only for either default, land or urban surface
           IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
-             CALL set_gas_flux( surf_def_h(0), emission_index_chem, chem_emission_att%units,       &
-                                dum_var_3d(:,:,in), time_factor(in) )
+             CALL set_gas_flux( surf_def_h(0), emission_index_chem, chem_emission_att%units,    &
+                                dum_var_3d, time_factor(in) )
           ELSE
-             CALL set_gas_flux( surf_usm_h, emission_index_chem, chem_emission_att%units,          &
-                                dum_var_3d(:,:,in), time_factor(in) )
-             CALL set_gas_flux( surf_lsm_h, emission_index_chem, chem_emission_att%units,          &
-                                dum_var_3d(:,:,in), time_factor(in) )
+             CALL set_gas_flux( surf_usm_h, emission_index_chem, chem_emission_att%units,       &
+                                dum_var_3d, time_factor(in) )
+             CALL set_gas_flux( surf_lsm_h, emission_index_chem, chem_emission_att%units,       &
+                                dum_var_3d, time_factor(in) )
           ENDIF
        ENDDO
@@ -9451,8 +9586,11 @@
     ELSEIF ( lod_gas_emissions == 2 )  THEN
 !
-!--    Obtain time index for current input starting at 0.
-!--    @todo: At the moment emission data and simulated time correspond to each other.
-       chem_emission_att%i_hour = MINLOC( ABS( gas_emission_time - time_since_reference_point ),   &
-                                          DIM = 1 ) - 1
+!--    Get time_utc_init from origin_date_time
+       CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
+!
+!--    Obtain time index for current point in time. Note, the time coordinate in the input file is
+!--    relative to time_utc_init.
+       chem_emission_att%i_hour = MINLOC( ABS( gas_emission_time - ( time_utc_init +               &
+                                         MAX( time_since_reference_point, 0.0_wp) ) ), DIM = 1 ) - 1
 !
 !--    Allocate the data input array always before reading in the data and deallocate after (NOTE 
@@ -9527,5 +9665,5 @@
        REAL(wp), DIMENSION(ngases_salsa) ::  conv     !< unit conversion factor
 
-       REAL(wp), DIMENSION(nys:nyn,nxl:nxr,chem_emission_att%n_emiss_species), INTENT(in) ::  source_array  !<
+       REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species), INTENT(in) ::  source_array  !<
 
        TYPE(surf_type), INTENT(inout) :: surface  !< respective surface type
@@ -11967,3 +12105,981 @@
  END FUNCTION
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set the lateral and top boundary conditions in case the PALM domain is
+!> nested offline in a mesoscale model. Further, average boundary data and
+!> determine mean profiles, further used for correct damping in the sponge
+!> layer.
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_nesting_offl_bc
+
+    USE control_parameters,                                                                        &
+        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, dt_3d,              &
+               time_since_reference_point
+
+    USE indices,                                                                                   &
+        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb, nzt
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i    !< running index x-direction
+    INTEGER(iwp) ::  ib   !< running index for aerosol number bins
+    INTEGER(iwp) ::  ic   !< running index for aerosol mass bins
+    INTEGER(iwp) ::  icc  !< running index for aerosol mass bins
+    INTEGER(iwp) ::  ig   !< running index for gaseous species
+    INTEGER(iwp) ::  j    !< running index y-direction
+    INTEGER(iwp) ::  k    !< running index z-direction
+
+    REAL(wp) ::  fac_dt  !< interpolation factor
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_mconc    !< reference profile for aerosol mass
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_mconc_l  !< reference profile for aerosol mass: subdomain
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_nconc    !< reference profile for aerosol number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_nconc_l  !< reference profile for aerosol_number: subdomain
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_gconc    !< reference profile for gases
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_gconc_l  !< reference profile for gases: subdomain
+
+!
+!-- Skip input if no forcing from larger-scale models is applied.
+    IF ( .NOT. nesting_offline_salsa )  RETURN
+!
+!-- Allocate temporary arrays to compute salsa mean profiles
+    ALLOCATE( ref_gconc(nzb:nzt+1,1:ngases_salsa), ref_gconc_l(nzb:nzt+1,1:ngases_salsa),          &
+              ref_mconc(nzb:nzt+1,1:nbins_aerosol*ncomponents_mass),                               &
+              ref_mconc_l(nzb:nzt+1,1:nbins_aerosol*ncomponents_mass),                             &
+              ref_nconc(nzb:nzt+1,1:nbins_aerosol), ref_nconc_l(nzb:nzt+1,1:nbins_aerosol) )
+    ref_gconc   = 0.0_wp
+    ref_gconc_l = 0.0_wp
+    ref_mconc   = 0.0_wp
+    ref_mconc_l = 0.0_wp
+    ref_nconc   = 0.0_wp
+    ref_nconc_l = 0.0_wp
+
+!
+!-- Determine interpolation factor and limit it to 1. This is because t+dt can slightly exceed
+!-- time(tind_p) before boundary data is updated again.
+    fac_dt = ( time_since_reference_point - salsa_nest_offl%time(salsa_nest_offl%tind)  + dt_3d )  &
+             / ( salsa_nest_offl%time(salsa_nest_offl%tind_p) -                                    &
+                 salsa_nest_offl%time(salsa_nest_offl%tind) )
+    fac_dt = MIN( 1.0_wp, fac_dt )
+
+    IF ( bc_dirichlet_l )  THEN
+       DO  ib = 1, nbins_aerosol
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                aerosol_number(ib)%conc(k,j,-1) = ( 1.0_wp - fac_dt ) *                            &
+                                                  salsa_nest_offl%nconc_left(0,k,j,ib) + fac_dt *  &
+                                                  salsa_nest_offl%nconc_left(1,k,j,ib)
+             ENDDO
+             ref_nconc_l(nzb+1:nzt,ib) = ref_nconc_l(nzb+1:nzt,ib) +                               &
+                                         aerosol_number(ib)%conc(nzb+1:nzt,j,-1)
+          ENDDO
+          DO  ic = 1, ncomponents_mass
+             icc = ( ic-1 ) * nbins_aerosol + ib
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   aerosol_mass(icc)%conc(k,j,-1) = ( 1.0_wp - fac_dt ) *                          &
+                                                    salsa_nest_offl%mconc_left(0,k,j,icc) + fac_dt &
+                                                    * salsa_nest_offl%mconc_left(1,k,j,icc)
+                ENDDO
+                ref_mconc_l(nzb+1:nzt,icc) = ref_mconc_l(nzb+1:nzt,icc) +                          &
+                                             aerosol_mass(icc)%conc(nzb+1:nzt,j,-1)
+             ENDDO
+          ENDDO
+       ENDDO
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          DO  ig = 1, ngases_salsa
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   salsa_gas(ig)%conc(k,j,-1) = ( 1.0_wp - fac_dt ) *                              &
+                                                salsa_nest_offl%gconc_left(0,k,j,ig) + fac_dt *    &
+                                                salsa_nest_offl%gconc_left(1,k,j,ig)
+                ENDDO
+                ref_gconc_l(nzb+1:nzt,ig) = ref_gconc_l(nzb+1:nzt,ig) +                            &
+                                            salsa_gas(ig)%conc(nzb+1:nzt,j,-1)
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+
+    IF ( bc_dirichlet_r )  THEN
+       DO  ib = 1, nbins_aerosol
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                aerosol_number(ib)%conc(k,j,nxr+1) = ( 1.0_wp - fac_dt ) *                         &
+                                                  salsa_nest_offl%nconc_right(0,k,j,ib) + fac_dt * &
+                                                  salsa_nest_offl%nconc_right(1,k,j,ib)
+             ENDDO
+             ref_nconc_l(nzb+1:nzt,ib) = ref_nconc_l(nzb+1:nzt,ib) +                               &
+                                         aerosol_number(ib)%conc(nzb+1:nzt,j,nxr+1)
+          ENDDO
+          DO  ic = 1, ncomponents_mass
+             icc = ( ic-1 ) * nbins_aerosol + ib
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   aerosol_mass(icc)%conc(k,j,nxr+1) = ( 1.0_wp - fac_dt ) *                       &
+                                                    salsa_nest_offl%mconc_right(0,k,j,icc) + fac_dt&
+                                                    * salsa_nest_offl%mconc_right(1,k,j,icc)
+                ENDDO
+                ref_mconc_l(nzb+1:nzt,icc) = ref_mconc_l(nzb+1:nzt,icc) +                          &
+                                             aerosol_mass(icc)%conc(nzb+1:nzt,j,nxr+1)
+             ENDDO
+          ENDDO
+       ENDDO
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          DO  ig = 1, ngases_salsa
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   salsa_gas(ig)%conc(k,j,nxr+1) = ( 1.0_wp - fac_dt ) *                           &
+                                                   salsa_nest_offl%gconc_right(0,k,j,ig) + fac_dt *&
+                                                   salsa_nest_offl%gconc_right(1,k,j,ig)
+                ENDDO
+                ref_gconc_l(nzb+1:nzt,ig) = ref_gconc_l(nzb+1:nzt,ig) +                            &
+                                            salsa_gas(ig)%conc(nzb+1:nzt,j,nxr+1)
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+
+    IF ( bc_dirichlet_n )  THEN
+       DO  ib = 1, nbins_aerosol
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                aerosol_number(ib)%conc(k,nyn+1,i) = ( 1.0_wp - fac_dt ) *                         &
+                                                  salsa_nest_offl%nconc_north(0,k,i,ib) + fac_dt * &
+                                                  salsa_nest_offl%nconc_north(1,k,i,ib)
+             ENDDO
+             ref_nconc_l(nzb+1:nzt,ib) = ref_nconc_l(nzb+1:nzt,ib) +                               &
+                                         aerosol_number(ib)%conc(nzb+1:nzt,nyn+1,i)
+          ENDDO
+          DO  ic = 1, ncomponents_mass
+             icc = ( ic-1 ) * nbins_aerosol + ib
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   aerosol_mass(icc)%conc(k,nyn+1,i) = ( 1.0_wp - fac_dt ) *                       &
+                                                    salsa_nest_offl%mconc_north(0,k,i,icc) + fac_dt&
+                                                    * salsa_nest_offl%mconc_north(1,k,i,icc)
+                ENDDO
+                ref_mconc_l(nzb+1:nzt,icc) = ref_mconc_l(nzb+1:nzt,icc) +                          &
+                                             aerosol_mass(icc)%conc(nzb+1:nzt,nyn+1,i)
+             ENDDO
+          ENDDO
+       ENDDO
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          DO  ig = 1, ngases_salsa
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   salsa_gas(ig)%conc(k,nyn+1,i) = ( 1.0_wp - fac_dt ) *                           &
+                                                   salsa_nest_offl%gconc_north(0,k,i,ig) + fac_dt *&
+                                                   salsa_nest_offl%gconc_north(1,k,i,ig)
+                ENDDO
+                ref_gconc_l(nzb+1:nzt,ig) = ref_gconc_l(nzb+1:nzt,ig) +                            &
+                                            salsa_gas(ig)%conc(nzb+1:nzt,nyn+1,i)
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+
+    IF ( bc_dirichlet_s )  THEN
+       DO  ib = 1, nbins_aerosol
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                aerosol_number(ib)%conc(k,-1,i) = ( 1.0_wp - fac_dt ) *                            &
+                                                  salsa_nest_offl%nconc_south(0,k,i,ib) + fac_dt * &
+                                                  salsa_nest_offl%nconc_south(1,k,i,ib)
+             ENDDO
+             ref_nconc_l(nzb+1:nzt,ib) = ref_nconc_l(nzb+1:nzt,ib) +                               &
+                                         aerosol_number(ib)%conc(nzb+1:nzt,-1,i)
+          ENDDO
+          DO  ic = 1, ncomponents_mass
+             icc = ( ic-1 ) * nbins_aerosol + ib
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   aerosol_mass(icc)%conc(k,-1,i) = ( 1.0_wp - fac_dt ) *                          &
+                                                    salsa_nest_offl%mconc_south(0,k,i,icc) + fac_dt&
+                                                    * salsa_nest_offl%mconc_south(1,k,i,icc)
+                ENDDO
+                ref_mconc_l(nzb+1:nzt,icc) = ref_mconc_l(nzb+1:nzt,icc) +                          &
+                                             aerosol_mass(icc)%conc(nzb+1:nzt,-1,i)
+             ENDDO
+          ENDDO
+       ENDDO
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          DO  ig = 1, ngases_salsa
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   salsa_gas(ig)%conc(k,-1,i) = ( 1.0_wp - fac_dt ) *                              &
+                                                salsa_nest_offl%gconc_south(0,k,i,ig) + fac_dt *   &
+                                                salsa_nest_offl%gconc_south(1,k,i,ig)
+                ENDDO
+                ref_gconc_l(nzb+1:nzt,ig) = ref_gconc_l(nzb+1:nzt,ig) +                            &
+                                            salsa_gas(ig)%conc(nzb+1:nzt,-1,i)
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+!
+!-- Top boundary
+    DO  ib = 1, nbins_aerosol
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             aerosol_number(ib)%conc(nzt+1,j,i) = ( 1.0_wp - fac_dt ) *                            &
+                                                  salsa_nest_offl%nconc_top(0,j,i,ib) + fac_dt *   &
+                                                  salsa_nest_offl%nconc_top(1,j,i,ib)
+             ref_nconc_l(nzt+1,ib) = ref_nconc_l(nzt+1,ib) + aerosol_number(ib)%conc(nzt+1,j,i)
+          ENDDO
+       ENDDO
+       DO  ic = 1, ncomponents_mass
+          icc = ( ic-1 ) * nbins_aerosol + ib
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                aerosol_mass(icc)%conc(nzt+1,j,i) = ( 1.0_wp - fac_dt ) *                          &
+                                                    salsa_nest_offl%mconc_top(0,j,i,icc) + fac_dt *&
+                                                    salsa_nest_offl%mconc_top(1,j,i,icc)
+                ref_mconc_l(nzt+1,icc) = ref_mconc_l(nzt+1,icc) + aerosol_mass(icc)%conc(nzt+1,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                salsa_gas(ig)%conc(nzt+1,j,i) = ( 1.0_wp - fac_dt ) *                              &
+                                                salsa_nest_offl%gconc_top(0,j,i,ig) + fac_dt *     &
+                                                salsa_nest_offl%gconc_top(1,j,i,ig)
+                ref_gconc_l(nzt+1,ig) = ref_gconc_l(nzt+1,ig) + salsa_gas(ig)%conc(nzt+1,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+!
+!-- Do local exchange
+    DO  ib = 1, nbins_aerosol
+       CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
+       DO  ic = 1, ncomponents_mass
+          icc = ( ic-1 ) * nbins_aerosol + ib
+          CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
+       ENDDO
+    ENDDO
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
+       ENDDO
+    ENDIF
+!
+!-- In case of Rayleigh damping, where the initial profiles are still used, update these profiles
+!-- from the averaged boundary data. But first, average these data.
+#if defined( __parallel )
+    IF ( .NOT. salsa_gases_from_chem )                                                             &
+       CALL MPI_ALLREDUCE( ref_gconc_l, ref_gconc, ( nzt+1-nzb+1 ) * SIZE( ref_gconc(nzb,:) ),     &
+                           MPI_REAL, MPI_SUM, comm2d, ierr )
+    CALL MPI_ALLREDUCE( ref_mconc_l, ref_mconc, ( nzt+1-nzb+1 ) * SIZE( ref_mconc(nzb,:) ),        &
+                        MPI_REAL, MPI_SUM, comm2d, ierr )
+    CALL MPI_ALLREDUCE( ref_nconc_l, ref_nconc, ( nzt+1-nzb+1 ) * SIZE( ref_nconc(nzb,:) ),        &
+                        MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+    IF ( .NOT. salsa_gases_from_chem )  ref_gconc = ref_gconc_l
+    ref_mconc = ref_mconc_l
+    ref_nconc = ref_nconc_l
+#endif
+!
+!-- Average data. Note, reference profiles up to nzt are derived from lateral boundaries, at the
+!-- model top it is derived from the top boundary. Thus, number of input data is different from
+!-- nzb:nzt compared to nzt+1.
+!-- Derived from lateral boundaries.
+    IF ( .NOT. salsa_gases_from_chem )                                                             &
+       ref_gconc(nzb:nzt,:) = ref_gconc(nzb:nzt,:) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), KIND = wp )
+    ref_mconc(nzb:nzt,:) = ref_mconc(nzb:nzt,:) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), KIND = wp )
+    ref_nconc(nzb:nzt,:) = ref_nconc(nzb:nzt,:) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), KIND = wp )
+!
+!-- Derived from top boundary
+    IF ( .NOT. salsa_gases_from_chem )                                                             &
+       ref_gconc(nzt+1,:) = ref_gconc(nzt+1,:) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
+    ref_mconc(nzt+1,:) = ref_mconc(nzt+1,:) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
+    ref_nconc(nzt+1,:) = ref_nconc(nzt+1,:) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
+!
+!-- Write onto init profiles, which are used for damping. Also set lower boundary condition.
+    DO  ib = 1, nbins_aerosol
+       aerosol_number(ib)%init(:)   = ref_nconc(:,ib)
+       aerosol_number(ib)%init(nzb) = aerosol_number(ib)%init(nzb+1)
+       DO  ic = 1, ncomponents_mass
+          icc = ( ic-1 ) * nbins_aerosol + ib
+          aerosol_mass(icc)%init(:)   = ref_mconc(:,icc)
+          aerosol_mass(icc)%init(nzb) = aerosol_mass(icc)%init(nzb+1)
+       ENDDO
+    ENDDO
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          salsa_gas(ig)%init(:)   = ref_gconc(:,ig)
+          salsa_gas(ig)%init(nzb) = salsa_gas(ig)%init(nzb+1)
+       ENDDO
+    ENDIF
+
+    DEALLOCATE( ref_gconc, ref_gconc_l, ref_mconc, ref_mconc_l, ref_nconc, ref_nconc_l )
+
+ END SUBROUTINE salsa_nesting_offl_bc
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Allocate arrays used to read boundary data from NetCDF file and initialize
+!> boundary data.
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_nesting_offl_init
+
+    USE control_parameters,                                                                        &
+        ONLY:  end_time, initializing_actions, spinup_time
+
+    USE palm_date_time_mod,                                                                        &
+        ONLY:  get_date_time
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  ib          !< running index for aerosol number bins
+    INTEGER(iwp) ::  ic          !< running index for aerosol mass bins
+    INTEGER(iwp) ::  icc         !< additional running index for aerosol mass bins
+    INTEGER(iwp) ::  ig          !< running index for gaseous species
+    INTEGER(iwp) ::  nmass_bins  !< number of aerosol mass bins
+
+    nmass_bins = nbins_aerosol * ncomponents_mass
+!
+!-- Get time_utc_init from origin_date_time
+    CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
+!
+!-- Allocate arrays for reading boundary values. Arrays will incorporate 2 time levels in order to
+!-- interpolate in between.
+    IF ( nesting_offline_salsa )  THEN
+       IF ( bc_dirichlet_l )  THEN
+          ALLOCATE( salsa_nest_offl%nconc_left(0:1,nzb+1:nzt,nys:nyn,1:nbins_aerosol) )
+          ALLOCATE( salsa_nest_offl%mconc_left(0:1,nzb+1:nzt,nys:nyn,1:nmass_bins) )
+       ENDIF
+       IF ( bc_dirichlet_r )  THEN
+          ALLOCATE( salsa_nest_offl%nconc_right(0:1,nzb+1:nzt,nys:nyn,1:nbins_aerosol) )
+          ALLOCATE( salsa_nest_offl%mconc_right(0:1,nzb+1:nzt,nys:nyn,1:nmass_bins) )
+       ENDIF
+       IF ( bc_dirichlet_n )  THEN
+          ALLOCATE( salsa_nest_offl%nconc_north(0:1,nzb+1:nzt,nxl:nxr,1:nbins_aerosol) )
+          ALLOCATE( salsa_nest_offl%mconc_north(0:1,nzb+1:nzt,nxl:nxr,1:nmass_bins) )
+       ENDIF
+       IF ( bc_dirichlet_s )  THEN
+          ALLOCATE( salsa_nest_offl%nconc_south(0:1,nzb+1:nzt,nxl:nxr,1:nbins_aerosol) )
+          ALLOCATE( salsa_nest_offl%mconc_south(0:1,nzb+1:nzt,nxl:nxr,1:nmass_bins) )
+       ENDIF
+       ALLOCATE( salsa_nest_offl%nconc_top(0:1,nys:nyn,nxl:nxr,1:nbins_aerosol) )
+       ALLOCATE( salsa_nest_offl%mconc_top(0:1,nys:nyn,nxl:nxr,1:nmass_bins) )
+
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          IF ( bc_dirichlet_l )  THEN
+             ALLOCATE( salsa_nest_offl%gconc_left(0:1,nzb+1:nzt,nys:nyn,1:ngases_salsa) )
+          ENDIF
+          IF ( bc_dirichlet_r )  THEN
+             ALLOCATE( salsa_nest_offl%gconc_right(0:1,nzb+1:nzt,nys:nyn,1:ngases_salsa) )
+          ENDIF
+          IF ( bc_dirichlet_n )  THEN
+             ALLOCATE( salsa_nest_offl%gconc_north(0:1,nzb+1:nzt,nxl:nxr,1:ngases_salsa) )
+          ENDIF
+          IF ( bc_dirichlet_s )  THEN
+             ALLOCATE( salsa_nest_offl%gconc_south(0:1,nzb+1:nzt,nxl:nxr,1:ngases_salsa) )
+          ENDIF
+          ALLOCATE( salsa_nest_offl%gconc_top(0:1,nys:nyn,nxl:nxr,1:ngases_salsa) )
+       ENDIF
+
+!
+!--    Read data at lateral and top boundaries from a larger-scale model
+       CALL salsa_nesting_offl_input
+!
+!--    Check if sufficient time steps are provided to cover the entire simulation. Note, dynamic
+!--    input is only required for the 3D simulation, not for the soil/wall spinup. However, as the
+!--    spinup time is added to the end_time, this must be considered here.
+       IF ( end_time - spinup_time >                                           &
+            salsa_nest_offl%time(salsa_nest_offl%nt-1) - time_utc_init )  THEN
+          message_string = 'end_time of the simulation exceeds the time dimension in the dynamic'//&
+                           ' input file.'
+          CALL message( 'salsa_nesting_offl_init', 'PA0681', 1, 2, 0, 6, 0 ) 
+       ENDIF
+
+       IF ( salsa_nest_offl%time(0) /= time_utc_init )  THEN
+          message_string = 'Offline nesting: time dimension must start at time_utc_init.'
+          CALL message( 'salsa_nesting_offl_init', 'PA0682', 1, 2, 0, 6, 0 )
+       ENDIF
+!
+!--    Initialize boundary data. Please note, do not initialize boundaries in case of restart runs.
+       IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  read_restart_data_salsa )  &
+       THEN
+          IF ( bc_dirichlet_l )  THEN
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc(nzb+1:nzt,nys:nyn,-1) =                                    &
+                                                 salsa_nest_offl%nconc_left(0,nzb+1:nzt,nys:nyn,ib)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc(nzb+1:nzt,nys:nyn,-1) =                                  &
+                                                 salsa_nest_offl%mconc_left(0,nzb+1:nzt,nys:nyn,icc)
+                ENDDO
+             ENDDO
+             DO  ig = 1, ngases_salsa
+                salsa_gas(ig)%conc(nzb+1:nzt,nys:nyn,-1) =                                         &
+                                                 salsa_nest_offl%gconc_left(0,nzb+1:nzt,nys:nyn,ig)
+             ENDDO
+          ENDIF
+          IF ( bc_dirichlet_r )  THEN
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc(nzb+1:nzt,nys:nyn,nxr+1) =                                 &
+                                                salsa_nest_offl%nconc_right(0,nzb+1:nzt,nys:nyn,ib)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc(nzb+1:nzt,nys:nyn,nxr+1) =                               &
+                                                salsa_nest_offl%mconc_right(0,nzb+1:nzt,nys:nyn,icc)
+                ENDDO
+             ENDDO
+             DO  ig = 1, ngases_salsa
+                salsa_gas(ig)%conc(nzb+1:nzt,nys:nyn,nxr+1) =                                      &
+                                                 salsa_nest_offl%gconc_right(0,nzb+1:nzt,nys:nyn,ig)
+             ENDDO
+          ENDIF
+          IF ( bc_dirichlet_n )  THEN
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =                                 &
+                                                salsa_nest_offl%nconc_north(0,nzb+1:nzt,nxl:nxr,ib)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =                               &
+                                                salsa_nest_offl%mconc_north(0,nzb+1:nzt,nxl:nxr,icc)
+                ENDDO
+             ENDDO
+             DO  ig = 1, ngases_salsa
+                salsa_gas(ig)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =                                      &
+                                                 salsa_nest_offl%gconc_north(0,nzb+1:nzt,nxl:nxr,ig)
+             ENDDO
+          ENDIF
+          IF ( bc_dirichlet_s )  THEN
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc(nzb+1:nzt,-1,nxl:nxr) =                                    &
+                                                salsa_nest_offl%nconc_south(0,nzb+1:nzt,nxl:nxr,ib)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc(nzb+1:nzt,-1,nxl:nxr) =                                  &
+                                                salsa_nest_offl%mconc_south(0,nzb+1:nzt,nxl:nxr,icc)
+                ENDDO
+             ENDDO
+             DO  ig = 1, ngases_salsa
+                salsa_gas(ig)%conc(nzb+1:nzt,-1,nxl:nxr) =                                         &
+                                                 salsa_nest_offl%gconc_south(0,nzb+1:nzt,nxl:nxr,ig)
+             ENDDO
+          ENDIF
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE salsa_nesting_offl_init
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set the lateral and top boundary conditions in case the PALM domain is
+!> nested offline in a mesoscale model. Further, average boundary data and
+!> determine mean profiles, further used for correct damping in the sponge
+!> layer.
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_nesting_offl_input
+
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
+               inquire_num_variables, inquire_variable_names,                                      &
+               get_dimension_length, open_read_file
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=25) ::  vname  !< variable name
+
+    INTEGER(iwp) ::  ic        !< running index for aerosol chemical components
+    INTEGER(iwp) ::  ig        !< running index for gases
+    INTEGER(iwp) ::  num_vars  !< number of variables in netcdf input file
+
+!
+!-- Skip input if no forcing from larger-scale models is applied.
+    IF ( .NOT. nesting_offline_salsa )  RETURN
+!
+!-- Initialise
+    IF ( .NOT. salsa_nest_offl%init )  THEN
+
+#if defined ( __netcdf )
+!
+!--    Open file in read-only mode
+       CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ),                   &
+                            salsa_nest_offl%id_dynamic )
+!
+!--    At first, inquire all variable names.
+       CALL inquire_num_variables( salsa_nest_offl%id_dynamic, num_vars )
+!
+!--    Allocate memory to store variable names.
+       ALLOCATE( salsa_nest_offl%var_names(1:num_vars) )
+       CALL inquire_variable_names( salsa_nest_offl%id_dynamic, salsa_nest_offl%var_names )
+!
+!--    Read time dimension, allocate memory and finally read time array
+       CALL get_dimension_length( salsa_nest_offl%id_dynamic, salsa_nest_offl%nt,&
+                                                    'time' )
+
+       IF ( check_existence( salsa_nest_offl%var_names, 'time' ) )  THEN
+          ALLOCATE( salsa_nest_offl%time(0:salsa_nest_offl%nt-1) )
+          CALL get_variable( salsa_nest_offl%id_dynamic, 'time', salsa_nest_offl%time )
+       ENDIF
+!
+!--    Read the vertical dimension
+       CALL get_dimension_length( salsa_nest_offl%id_dynamic,                    &
+                                                    salsa_nest_offl%nzu, 'z' )
+       ALLOCATE( salsa_nest_offl%zu_atmos(1:salsa_nest_offl%nzu) )
+       CALL get_variable( salsa_nest_offl%id_dynamic, 'z', salsa_nest_offl%zu_atmos )
+!
+!--    Read the number of aerosol chemical components
+       CALL get_dimension_length( salsa_nest_offl%id_dynamic,                    &
+                                                    salsa_nest_offl%ncc, 'composition_index' )
+!
+!--    Read the names of aerosol chemical components
+       CALL get_variable( salsa_nest_offl%id_dynamic, 'composition_name', salsa_nest_offl%cc_name, &
+                          salsa_nest_offl%ncc )
+!
+!--    Define the index of each chemical component in the model
+       DO  ic = 1, salsa_nest_offl%ncc
+          SELECT CASE ( TRIM( salsa_nest_offl%cc_name(ic) ) )
+             CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' )
+                salsa_nest_offl%cc_in2mod(1) = ic
+             CASE ( 'OC', 'oc' )
+                salsa_nest_offl%cc_in2mod(2) = ic
+             CASE ( 'BC', 'bc' )
+                salsa_nest_offl%cc_in2mod(3) = ic
+             CASE ( 'DU', 'du' )
+                salsa_nest_offl%cc_in2mod(4) = ic
+             CASE ( 'SS', 'ss' )
+                salsa_nest_offl%cc_in2mod(5) = ic
+             CASE ( 'HNO3', 'hno3', 'NO3', 'no3', 'NO', 'no' )
+                salsa_nest_offl%cc_in2mod(6) = ic
+             CASE ( 'NH3', 'nh3', 'NH4', 'nh4', 'NH', 'nh' )
+                salsa_nest_offl%cc_in2mod(7) = ic
+          END SELECT
+       ENDDO
+       IF ( SUM( salsa_nest_offl%cc_in2mod ) == 0 )  THEN
+          message_string = 'None of the aerosol chemical components in ' //                        &
+                           TRIM( input_file_dynamic ) // ' correspond to ones applied in SALSA.'
+          CALL message( 'salsa_mod: salsa_nesting_offl_input',                      &
+                        'PA0662', 2, 2, 0, 6, 0 )
+       ENDIF
+#endif
+    ENDIF
+!
+!-- Check if dynamic driver data input is required.
+    IF ( salsa_nest_offl%time(salsa_nest_offl%tind_p) <= MAX( time_since_reference_point, 0.0_wp)  &
+         + time_utc_init  .OR.  .NOT.  salsa_nest_offl%init )  THEN
+       CONTINUE
+!
+!-- Return otherwise
+    ELSE
+       RETURN
+    ENDIF
+!
+!-- Obtain time index for current point in time.
+    salsa_nest_offl%tind = MINLOC( ABS( salsa_nest_offl%time - ( time_utc_init +                   &
+                                        MAX( time_since_reference_point, 0.0_wp) ) ), DIM = 1 ) - 1
+    salsa_nest_offl%tind_p = salsa_nest_offl%tind + 1
+!
+!-- Open file in read-only mode
+#if defined ( __netcdf )
+
+    CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ),                      &
+                         salsa_nest_offl%id_dynamic )
+!
+!-- Read data at the western boundary
+    CALL get_variable( salsa_nest_offl%id_dynamic, 'ls_forcing_left_aerosol',                      &
+                       salsa_nest_offl%nconc_left,                                                 &
+                       MERGE( 0, 1, bc_dirichlet_l ), MERGE( nbins_aerosol-1, 0, bc_dirichlet_l ), &
+                       MERGE( nys, 1, bc_dirichlet_l ), MERGE( nyn, 0, bc_dirichlet_l ),           &
+                       MERGE( nzb, 1, bc_dirichlet_l ), MERGE( nzt-1, 0, bc_dirichlet_l ),         &
+                       MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_l ),                         &
+                       MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_l  ) )
+    IF ( bc_dirichlet_l )  THEN
+       salsa_nest_offl%nconc_left = MAX( nclim, salsa_nest_offl%nconc_left )
+       CALL nesting_offl_aero_mass( salsa_nest_offl%tind, salsa_nest_offl%tind_p, nzb+1, nzt, nys, &
+                                    nyn, 'ls_forcing_left_mass_fracs_a', 1 )
+    ENDIF
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          vname = salsa_nest_offl%char_l // salsa_nest_offl%gas_name(ig)
+          CALL get_variable( salsa_nest_offl%id_dynamic, TRIM( vname ),                            &
+                             salsa_nest_offl%gconc_left(:,:,:,ig),                                 &
+                             MERGE( nys, 1, bc_dirichlet_l ), MERGE( nyn, 0, bc_dirichlet_l ),     &
+                             MERGE( nzb, 1, bc_dirichlet_l ), MERGE( nzt-1, 0, bc_dirichlet_l ),   &
+                             MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_l ),                   &
+                             MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_l ) )
+          IF ( bc_dirichlet_l )  salsa_nest_offl%gconc_left(:,:,:,ig) =                            &
+                                                  MAX( nclim, salsa_nest_offl%gconc_left(:,:,:,ig) )
+       ENDDO
+    ENDIF
+!
+!-- Read data at the eastern boundary
+    CALL get_variable( salsa_nest_offl%id_dynamic, 'ls_forcing_right_aerosol',                     &
+                       salsa_nest_offl%nconc_right,                                                &
+                       MERGE( 0, 1, bc_dirichlet_r ), MERGE( nbins_aerosol-1, 0, bc_dirichlet_r ), &
+                       MERGE( nys, 1, bc_dirichlet_r ), MERGE( nyn, 0, bc_dirichlet_r ),           &
+                       MERGE( nzb, 1, bc_dirichlet_r ), MERGE( nzt-1, 0, bc_dirichlet_r ),         &
+                       MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_r ),                         &
+                       MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_r ) )
+    IF ( bc_dirichlet_r )  THEN
+       salsa_nest_offl%nconc_right = MAX( nclim, salsa_nest_offl%nconc_right )
+       CALL nesting_offl_aero_mass( salsa_nest_offl%tind, salsa_nest_offl%tind_p, nzb+1, nzt, nys, &
+                                    nyn, 'ls_forcing_right_mass_fracs_a', 2 )
+    ENDIF
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          vname = salsa_nest_offl%char_r // salsa_nest_offl%gas_name(ig)
+          CALL get_variable( salsa_nest_offl%id_dynamic, TRIM( vname ),                            &
+                             salsa_nest_offl%gconc_right(:,:,:,ig),                                &
+                             MERGE( nys, 1, bc_dirichlet_r ), MERGE( nyn, 0, bc_dirichlet_r ),     &
+                             MERGE( nzb, 1, bc_dirichlet_r ), MERGE( nzt-1, 0, bc_dirichlet_r ),   &
+                             MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_r ),                   &
+                             MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_r ) )
+          IF ( bc_dirichlet_r )  salsa_nest_offl%gconc_right(:,:,:,ig) =                           &
+                                                 MAX( nclim, salsa_nest_offl%gconc_right(:,:,:,ig) )
+       ENDDO
+    ENDIF
+!
+!-- Read data at the northern boundary
+    CALL get_variable( salsa_nest_offl%id_dynamic, 'ls_forcing_north_aerosol',                     &
+                       salsa_nest_offl%nconc_north,                                                &
+                       MERGE( 0, 1, bc_dirichlet_n ), MERGE( nbins_aerosol-1, 0, bc_dirichlet_n ), &
+                       MERGE( nxl, 1, bc_dirichlet_n ), MERGE( nxr, 0, bc_dirichlet_n ),           &
+                       MERGE( nzb, 1, bc_dirichlet_n ), MERGE( nzt-1, 0, bc_dirichlet_n ),         &
+                       MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_n ),                         &
+                       MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_n ) )
+    IF ( bc_dirichlet_n )  THEN
+       salsa_nest_offl%nconc_north = MAX( nclim, salsa_nest_offl%nconc_north )
+       CALL nesting_offl_aero_mass( salsa_nest_offl%tind, salsa_nest_offl%tind_p, nzb+1, nzt, nxl, &
+                                    nxr, 'ls_forcing_north_mass_fracs_a', 3 )
+    ENDIF
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          vname = salsa_nest_offl%char_n // salsa_nest_offl%gas_name(ig)
+          CALL get_variable( salsa_nest_offl%id_dynamic, TRIM( vname ),                            &
+                             salsa_nest_offl%gconc_north(:,:,:,ig),                                &
+                             MERGE( nxl, 1, bc_dirichlet_n ), MERGE( nxr, 0, bc_dirichlet_n ),     &
+                             MERGE( nzb, 1, bc_dirichlet_n ), MERGE( nzt-1, 0, bc_dirichlet_n ),   &
+                             MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_n ),                   &
+                             MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_n ) )
+          IF ( bc_dirichlet_n )  salsa_nest_offl%gconc_north(:,:,:,ig) =                           &
+                                                 MAX( nclim, salsa_nest_offl%gconc_north(:,:,:,ig) )
+       ENDDO
+    ENDIF
+!
+!-- Read data at the southern boundary
+    CALL get_variable( salsa_nest_offl%id_dynamic, 'ls_forcing_south_aerosol',                     &
+                       salsa_nest_offl%nconc_south,                                                &
+                       MERGE( 0, 1, bc_dirichlet_s ), MERGE( nbins_aerosol-1, 0, bc_dirichlet_s ), &
+                       MERGE( nxl, 1, bc_dirichlet_s ), MERGE( nxr, 0, bc_dirichlet_s ),           &
+                       MERGE( nzb, 1, bc_dirichlet_s ), MERGE( nzt-1, 0, bc_dirichlet_s ),         &
+                       MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_s ),                         &
+                       MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_s ) )
+    IF ( bc_dirichlet_s )  THEN
+       salsa_nest_offl%nconc_south = MAX( nclim, salsa_nest_offl%nconc_south )
+       CALL nesting_offl_aero_mass( salsa_nest_offl%tind, salsa_nest_offl%tind_p, nzb+1, nzt, nxl, &
+                                    nxr, 'ls_forcing_south_mass_fracs_a', 4 )
+    ENDIF
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          vname = salsa_nest_offl%char_s // salsa_nest_offl%gas_name(ig)
+          CALL get_variable( salsa_nest_offl%id_dynamic, TRIM( vname ),                            &
+                             salsa_nest_offl%gconc_south(:,:,:,ig),                                &
+                             MERGE( nxl, 1, bc_dirichlet_s ), MERGE( nxr, 0, bc_dirichlet_s ),     &
+                             MERGE( nzb, 1, bc_dirichlet_s ), MERGE( nzt-1, 0, bc_dirichlet_s ),   &
+                             MERGE( salsa_nest_offl%tind,   1, bc_dirichlet_s ),                   &
+                             MERGE( salsa_nest_offl%tind_p, 0, bc_dirichlet_s ) )
+          IF ( bc_dirichlet_s )  salsa_nest_offl%gconc_south(:,:,:,ig) =                           &
+                                                 MAX( nclim, salsa_nest_offl%gconc_south(:,:,:,ig) )
+       ENDDO
+    ENDIF
+!
+!-- Read data at the top boundary
+    CALL get_variable( salsa_nest_offl%id_dynamic, 'ls_forcing_top_aerosol',                       &
+                       salsa_nest_offl%nconc_top(0:1,nys:nyn,nxl:nxr,1:nbins_aerosol),             &
+                       0, nbins_aerosol-1, nxl, nxr, nys, nyn, salsa_nest_offl%tind,               &
+                       salsa_nest_offl%tind_p )
+    salsa_nest_offl%nconc_top = MAX( nclim, salsa_nest_offl%nconc_top )
+    CALL nesting_offl_aero_mass( salsa_nest_offl%tind, salsa_nest_offl%tind_p, nys, nyn, nxl, nxr, &
+                                 'ls_forcing_top_mass_fracs_a', 5 )
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  ig = 1, ngases_salsa
+          vname = salsa_nest_offl%char_t // salsa_nest_offl%gas_name(ig)
+          CALL get_variable( salsa_nest_offl%id_dynamic, TRIM( vname ),                            &
+                             salsa_nest_offl%gconc_top(:,:,:,ig), nxl, nxr, nys, nyn,              &
+                             salsa_nest_offl%tind, salsa_nest_offl%tind_p )
+          salsa_nest_offl%gconc_top(:,:,:,ig) = MAX( nclim, salsa_nest_offl%gconc_top(:,:,:,ig) )
+       ENDDO
+    ENDIF
+!
+!-- Close input file
+    CALL close_input_file( salsa_nest_offl%id_dynamic )
+
+#endif
+!
+!-- Set control flag to indicate that initialization is already done
+    salsa_nest_offl%init = .TRUE.
+
+ END SUBROUTINE salsa_nesting_offl_input
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sets the mass concentrations to aerosol arrays in 2a and 2b.
+!------------------------------------------------------------------------------!
+ SUBROUTINE nesting_offl_aero_mass( ts, te, ks, ke, is, ie, varname_a, ibound )
+
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  get_variable
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=25) ::  varname_b  !< name for bins b
+
+    CHARACTER(LEN=*), INTENT(in) ::  varname_a  !< name for bins a
+
+    INTEGER(iwp) ::  ee                !< loop index: end
+    INTEGER(iwp) ::  i                 !< loop index
+    INTEGER(iwp) ::  ib                !< loop index
+    INTEGER(iwp) ::  ic                !< loop index
+    INTEGER(iwp) ::  k                 !< loop index
+    INTEGER(iwp) ::  ss                !< loop index: start
+    INTEGER(iwp) ::  t                 !< loop index
+    INTEGER(iwp) ::  type_so4_oc = -1  !<
+
+    INTEGER(iwp), INTENT(in) ::  ibound  !< index: 1=left, 2=right, 3=north, 4=south, 5=top
+    INTEGER(iwp), INTENT(in) ::  ie      !< loop index
+    INTEGER(iwp), INTENT(in) ::  is      !< loop index
+    INTEGER(iwp), INTENT(in) ::  ks      !< loop index
+    INTEGER(iwp), INTENT(in) ::  ke      !< loop index
+    INTEGER(iwp), INTENT(in) ::  ts      !< loop index
+    INTEGER(iwp), INTENT(in) ::  te      !< loop index
+
+    INTEGER(iwp), DIMENSION(maxspec) ::  cc_i2m   !<
+
+    REAL(wp) ::  pmf1a !< mass fraction in 1a
+
+    REAL(wp), DIMENSION(nbins_aerosol) ::  core   !< size of the bin mid aerosol particle
+
+    REAL(wp), DIMENSION(0:1,ks:ke,is:ie,1:nbins_aerosol) ::  to_nconc                   !<
+    REAL(wp), DIMENSION(0:1,ks:ke,is:ie,1:nbins_aerosol*ncomponents_mass) ::  to_mconc  !<
+
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mf2a !< Mass distributions for a
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  mf2b !< and b bins
+
+!
+!-- Variable name for insoluble mass fraction
+    varname_b = varname_a(1:LEN( TRIM( varname_a ) ) - 1 ) // 'b'
+!
+!-- Bin mean aerosol particle volume (m3)
+    core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3
+!
+!-- Allocate and read mass fraction arrays
+    ALLOCATE( mf2a(0:1,ks:ke,is:ie,1:salsa_nest_offl%ncc),                                         &
+              mf2b(0:1,ks:ke,is:ie,1:salsa_nest_offl%ncc) )
+    IF ( ibound == 5 )  THEN
+       CALL get_variable( salsa_nest_offl%id_dynamic, varname_a,                                   &
+                          mf2a(0:1,ks:ke,is:ie,1:salsa_nest_offl%ncc), 0, salsa_nest_offl%ncc-1,   &
+                          is, ie, ks, ke, ts, te )
+    ELSE
+       CALL get_variable( salsa_nest_offl%id_dynamic, varname_a,                                   &
+                          mf2a(0:1,ks:ke,is:ie,1:salsa_nest_offl%ncc), 0, salsa_nest_offl%ncc-1,   &
+                          is, ie, ks-1, ke-1, ts, te )
+    ENDIF
+!
+!-- If the chemical component is not activated, set its mass fraction to 0 to avoid mass inbalance
+    cc_i2m = salsa_nest_offl%cc_in2mod
+    IF ( index_so4 < 0  .AND. cc_i2m(1) > 0 )  mf2a(:,:,:,cc_i2m(1)) = 0.0_wp
+    IF ( index_oc < 0   .AND. cc_i2m(2) > 0 )  mf2a(:,:,:,cc_i2m(2)) = 0.0_wp
+    IF ( index_bc < 0   .AND. cc_i2m(3) > 0 )  mf2a(:,:,:,cc_i2m(3)) = 0.0_wp
+    IF ( index_du < 0   .AND. cc_i2m(4) > 0 )  mf2a(:,:,:,cc_i2m(4)) = 0.0_wp
+    IF ( index_ss < 0   .AND. cc_i2m(5) > 0 )  mf2a(:,:,:,cc_i2m(5)) = 0.0_wp
+    IF ( index_no < 0   .AND. cc_i2m(6) > 0 )  mf2a(:,:,:,cc_i2m(6)) = 0.0_wp
+    IF ( index_nh < 0   .AND. cc_i2m(7) > 0 )  mf2a(:,:,:,cc_i2m(7)) = 0.0_wp
+    mf2b = 0.0_wp
+!
+!-- Initialise variable type_so4_oc to indicate whether SO4 and/OC is included in mass fraction data
+    IF ( ( cc_i2m(1) > 0  .AND.  index_so4 > 0 )  .AND. ( cc_i2m(2) > 0  .AND.  index_oc > 0 ) )   &
+    THEN
+       type_so4_oc = 1
+    ELSEIF ( cc_i2m(1) > 0  .AND.  index_so4 > 0 )  THEN
+       type_so4_oc = 2
+    ELSEIF ( cc_i2m(2) > 0  .AND.  index_oc > 0 )  THEN
+       type_so4_oc = 3
+    ENDIF
+
+    SELECT CASE ( ibound )
+       CASE( 1 )
+          to_nconc = salsa_nest_offl%nconc_left
+          to_mconc = salsa_nest_offl%mconc_left
+       CASE( 2 )
+          to_nconc = salsa_nest_offl%nconc_right
+          to_mconc = salsa_nest_offl%mconc_right
+       CASE( 3 )
+          to_nconc = salsa_nest_offl%nconc_north
+          to_mconc = salsa_nest_offl%mconc_north
+       CASE( 4 )
+          to_nconc = salsa_nest_offl%nconc_south
+          to_mconc = salsa_nest_offl%mconc_south
+       CASE( 5 )
+          to_nconc = salsa_nest_offl%nconc_top
+          to_mconc = salsa_nest_offl%mconc_top
+    END SELECT
+!
+!-- Set mass concentrations:
+!
+!-- Regime 1:
+    SELECT CASE ( type_so4_oc )
+       CASE ( 1 )  ! Both SO4 and OC given
+
+          ss = ( index_so4 - 1 ) * nbins_aerosol + start_subrange_1a  ! start
+          ee = ( index_so4 - 1 ) * nbins_aerosol + end_subrange_1a    ! end
+          ib = start_subrange_1a
+          DO  ic = ss, ee
+             DO i = is, ie
+                DO k = ks, ke
+                   DO t = 0, 1
+                      pmf1a = mf2a(t,k,i,cc_i2m(1)) / ( mf2a(t,k,i,cc_i2m(1)) + mf2a(t,k,i,cc_i2m(2)) )
+                      to_mconc(t,k,i,ic) = pmf1a * to_nconc(t,k,i,ib) * core(ib) * arhoh2so4
+                   ENDDO
+                ENDDO
+             ENDDO
+             ib = ib + 1
+          ENDDO
+          ss = ( index_oc - 1 ) * nbins_aerosol + start_subrange_1a ! start
+          ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a   ! end
+          ib = start_subrange_1a
+          DO  ic = ss, ee
+             DO i = is, ie
+                DO k = ks, ke
+                   DO t = 0, 1
+                      pmf1a = mf2a(t,k,i,cc_i2m(2)) / ( mf2a(t,k,i,cc_i2m(1)) + mf2a(t,k,i,cc_i2m(2)) )
+                      to_mconc(t,k,i,ic) = pmf1a * to_nconc(t,k,i,ib) * core(ib) * arhooc
+                   ENDDO
+                ENDDO
+             ENDDO
+             ib = ib + 1
+          ENDDO
+       CASE ( 2 )  ! Only SO4
+          ss = ( index_so4 - 1 ) * nbins_aerosol + start_subrange_1a  ! start
+          ee = ( index_so4 - 1 ) * nbins_aerosol + end_subrange_1a    ! end
+          ib = start_subrange_1a
+          DO  ic = ss, ee
+             DO i = is, ie
+                DO k = ks, ke
+                   DO t = 0, 1
+                      to_mconc(t,k,i,ic) = to_nconc(t,k,i,ib) * core(ib) * arhoh2so4
+                   ENDDO
+                ENDDO
+             ENDDO
+             ib = ib + 1
+          ENDDO
+       CASE ( 3 )  ! Only OC
+          ss = ( index_oc - 1 ) * nbins_aerosol + start_subrange_1a ! start
+          ee = ( index_oc - 1 ) * nbins_aerosol + end_subrange_1a   ! end
+          ib = start_subrange_1a
+          DO  ic = ss, ee
+             DO i = is, ie
+                DO k = ks, ke
+                   DO t = 0, 1
+                      to_mconc(t,k,i,ic) = to_nconc(t,k,i,ib) * core(ib) * arhooc
+                   ENDDO
+                ENDDO
+             ENDDO
+             ib = ib + 1
+          ENDDO
+    END SELECT
+!
+!-- Regimes 2a and 2b:
+    IF ( index_so4 > 0 ) THEN
+       CALL set_nest_mass( index_so4, 1, arhoh2so4 )
+    ENDIF
+    IF ( index_oc > 0 ) THEN
+       CALL set_nest_mass( index_oc, 2, arhooc )
+    ENDIF
+    IF ( index_bc > 0 ) THEN
+       CALL set_nest_mass( index_bc, 3, arhobc )
+    ENDIF
+    IF ( index_du > 0 ) THEN
+       CALL set_nest_mass( index_du, 4, arhodu )
+    ENDIF
+    IF ( index_ss > 0 ) THEN
+       CALL set_nest_mass( index_ss, 5, arhoss )
+    ENDIF
+    IF ( index_no > 0 ) THEN
+       CALL set_nest_mass( index_no, 6, arhohno3 )
+    ENDIF
+    IF ( index_nh > 0 ) THEN
+       CALL set_nest_mass( index_nh, 7, arhonh3 )
+    ENDIF
+
+    DEALLOCATE( mf2a, mf2b )
+
+    SELECT CASE ( ibound )
+       CASE( 1 )
+          salsa_nest_offl%mconc_left = to_mconc
+       CASE( 2 )
+          salsa_nest_offl%mconc_right = to_mconc
+       CASE( 3 )
+          salsa_nest_offl%mconc_north = to_mconc
+       CASE( 4 )
+          salsa_nest_offl%mconc_south = to_mconc
+       CASE( 5 )
+          salsa_nest_offl%mconc_top = to_mconc
+    END SELECT
+
+    CONTAINS
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set nesting boundaries for aerosol mass.
+!------------------------------------------------------------------------------!
+    SUBROUTINE set_nest_mass( ispec, ispec_def, prho )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  ic   !< chemical component index: default
+       INTEGER(iwp) ::  icc  !< loop index: mass bin
+
+       INTEGER(iwp), INTENT(in) ::  ispec      !< aerosol species index
+       INTEGER(iwp), INTENT(in) ::  ispec_def  !< default aerosol species index
+
+       REAL(wp), INTENT(in) ::  prho !< aerosol density
+!
+!--    Define the index of the chemical component in the input data
+       ic = salsa_nest_offl%cc_in2mod(ispec_def)
+
+       DO i = is, ie
+          DO k = ks, ke
+             DO t = 0, 1
+!
+!--             Regime 2a:
+                ss = ( ispec - 1 ) * nbins_aerosol + start_subrange_2a
+                ee = ( ispec - 1 ) * nbins_aerosol + end_subrange_2a
+                ib = start_subrange_2a
+                DO icc = ss, ee
+                   to_mconc(t,k,i,icc) = MAX( 0.0_wp, mf2a(t,k,i,ic) / SUM( mf2a(t,k,i,:) ) ) *    &
+                                         to_nconc(t,k,i,ib) * core(ib) * prho
+                   ib = ib + 1
+                ENDDO
+!
+!--             Regime 2b:
+                IF ( .NOT. no_insoluble )  THEN
+!
+!--                 TODO!
+                    mf2b(t,k,i,ic) = mf2b(t,k,i,ic)
+                ENDIF
+             ENDDO   ! k
+
+          ENDDO   ! j
+       ENDDO   ! i
+
+    END SUBROUTINE set_nest_mass
+
+ END SUBROUTINE nesting_offl_aero_mass
+
+
  END MODULE salsa_mod
