Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 3579)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -317,5 +318,5 @@
        USE control_parameters,                                                 &
            ONLY:  humidity, loop_optimization, passive_scalar, ocean_mode,     &
-                  rans_tke_e, ws_scheme_mom, ws_scheme_sca
+                  rans_tke_e, ws_scheme_mom, ws_scheme_sca, salsa
 
        USE indices,                                                            &
@@ -399,8 +400,9 @@
              sums_wssas_ws_l = 0.0_wp
           ENDIF        
-!          
-!--       Mona: Now always allocated
-          ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
-          sums_salsa_ws_l = 0.0_wp         
+          
+          IF ( salsa )  THEN
+             ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_salsa_ws_l = 0.0_wp   
+          ENDIF
 
        ENDIF
@@ -1115,5 +1117,5 @@
        USE control_parameters,                                                 &
            ONLY:  humidity, passive_scalar, ocean_mode, ws_scheme_mom,         &
-                  ws_scheme_sca
+                  ws_scheme_sca, salsa
 
        USE kinds  
@@ -1158,8 +1160,6 @@
           ENDIF
           IF ( ocean_mode )  sums_wssas_ws_l = 0.0_wp
-       
-          sums_salsa_ws_l = 0.0_wp
-         
-          
+          IF ( salsa )       sums_salsa_ws_l = 0.0_wp
+
        ENDIF
 
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 3579)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 3582)
@@ -20,6 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
+!
 ! Former revisions:
 ! -----------------
@@ -177,5 +178,5 @@
     USE control_parameters,                                                    &
         ONLY:  air_chemistry, average_count_3d, biometeorology, doav, doav_n,  &
-               land_surface, ocean_mode, urban_surface, varnamelength
+               land_surface, ocean_mode, salsa, urban_surface, varnamelength
 
     USE cpulog,                                                                &
@@ -200,5 +201,5 @@
            
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_3d_data_averaging 
+        ONLY:  salsa_3d_data_averaging 
 
     USE turbulence_closure_mod,                                                &
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 3579)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -225,9 +226,9 @@
                bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
                bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val,            &
-               child_domain, constant_diffusion, coupling_mode, &
-               dt_3d, humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, &
-               ibc_s_t, ibc_uv_b, ibc_uv_t,                           &
-               intermediate_timestep_count, nesting_offline, nudging,          &
-               ocean_mode, passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax
+               child_domain, constant_diffusion, coupling_mode, dt_3d,         &
+               humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b,        &
+               ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count,       &
+               nesting_offline, nudging, ocean_mode, passive_scalar, rans_mode,&
+               rans_tke_e, tsc, salsa, use_cmax
 
     USE grid_variables,                                                        &
@@ -248,5 +249,5 @@
  
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_boundary_conds        
+        ONLY:  salsa_boundary_conds        
 
     USE surface_mod,                                                           &
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 3579)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 3582)
@@ -20,4 +20,6 @@
 ! Current revisions:
 ! -----------------
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! 
@@ -813,5 +815,5 @@
                             
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_check_data_output, salsa_check_parameters        
+        ONLY:  salsa_check_data_output, salsa_check_parameters        
 
     USE spectra_mod,                                                           &
@@ -821,6 +823,4 @@
 
     USE subsidence_mod
-
-    USE statistics
 
     USE surface_output_mod,                                                    &
Index: palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3579)
+++ palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3582)
@@ -22,5 +22,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! - Removed salsa dependency.
+! - Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but 
+!   salsa is (M. Kurppa)
 ! 
 ! Former revisions:
@@ -890,6 +892,4 @@
    USE indices,                                                               &
        ONLY: nnx,nny,nnz 
-   USE salsa_mod,                                                             &
-       ONLY:  salsa
    USE surface_mod,                                                           &
        ONLY:  surf_lsm_h,surf_def_h,surf_usm_h
@@ -1462,74 +1462,4 @@
 
                 ENDDO
-             ELSEIF ( salsa ) THEN
-                DO  m = 1, surf_def_h(0)%ns
-                   i = surf_def_h(0)%i(m)
-                   j = surf_def_h(0)%j(m)
-                   k = surf_def_h(0)%k(m)
-
-
-                   IF ( street_type_f%var(j,i) >= main_street_id  .AND.        &
-                        street_type_f%var(j,i) < max_street_id )               &
-                   THEN
-
-                      !> Cycle over already matched species
-                      DO  ispec=1,nspec_out
-
-                         !> PMs are already in mass units:micrograms: have to be converted to kilograms 
-                         IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1"       &
-                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
-                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")&
-                         THEN
-                                  
-                            surf_def_h(0)%cssws(match_spec_model(ispec),m) =   &
-                                  emiss_factor_main(match_spec_input(ispec)) * &
-                                  emis_distribution(1,j,i,ispec) * rho_air(k) /&
-                                  time_factor(1)
-                         ELSE
-
-                         !> Other Species: inputs are micromoles: have to be converted              
-                            surf_def_h(0)%cssws(match_spec_model(ispec),m) =   &
-                               emiss_factor_main(match_spec_input(ispec)) *    &
-                               emis_distribution(1,j,i,ispec) *                &
-                               conv_to_ratio(k,j,i) *  rho_air(k) / time_factor(1)
-                         ENDIF
-                      ENDDO
-
-                   ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.    &
-                            street_type_f%var(j,i) < main_street_id )          &
-                   THEN
-
-                   !> Cycle over already matched species
-                      DO  ispec=1,nspec_out
-
-                         !> PMs are already in mass units: micrograms
-                         IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1"   &
-                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" &
-                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10")&
-                         THEN
-
-                            surf_def_h(0)%cssws(match_spec_model(ispec),m) =   &
-                               emiss_factor_side(match_spec_input(ispec)) *    &
-                               emis_distribution(1,j,i,ispec) * rho_air(k) /   &
-                               time_factor(1)  
-                         ELSE
-               
-                            surf_def_h(0)%cssws(match_spec_model(ispec),m) =   &
-                               emiss_factor_side(match_spec_input(ispec)) *    &
-                               emis_distribution(1,j,i,ispec) *                &
-                               conv_to_ratio(k,j,i) * rho_air(k) / time_factor(1)
-                         ENDIF
-
-                      ENDDO
-
-                   ELSE
-
-                   !> If no street type is defined, then assign null emissions to all the species
-                      surf_def_h(0)%cssws(:,m) = 0.0_wp
-
-                   ENDIF
-
-                ENDDO
-
              ENDIF
 
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 3579)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -294,5 +295,5 @@
                ibc_uv_b, io_blocks, io_group, land_surface, message_string,    &
                ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
-               ocean_mode, psolver, section, simulated_time,                   &
+               ocean_mode, psolver, salsa, section, simulated_time,            &
                time_since_reference_point
 
@@ -334,5 +335,5 @@
    
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_data_output_2d    
+        ONLY:  salsa_data_output_2d    
 
     USE surface_mod,                                                           &
@@ -1353,6 +1354,7 @@
                 
                 IF ( .NOT. found  .AND.  salsa )  THEN
-                   CALL salsa_data_output_2d( av, do2d(av,ivar), found, grid,    &
-                                              mode, local_pf, two_d )
+                   CALL salsa_data_output_2d( av, do2d(av,ivar), found, grid,  &
+                                              mode, local_pf, two_d, nzb_do,   &
+                                              nzt_do)
                 ENDIF                 
 
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 3579)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 3582)
@@ -20,6 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
+!
 ! Former revisions:
 ! -----------------
@@ -254,5 +255,5 @@
                io_blocks, io_group, land_surface, message_string,              &
                ntdim_3d, nz_do3d,  ocean_mode, plant_canopy,                   &
-               psolver, simulated_time, time_since_reference_point,            &
+               psolver, salsa, simulated_time, time_since_reference_point,     &
                urban_surface, varnamelength
 
@@ -302,5 +303,5 @@
                
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_data_output_3d     
+        ONLY:  salsa_data_output_3d     
 
     USE turbulence_closure_mod,                                                &
@@ -791,5 +792,6 @@
 !--          SALSA output
              IF ( .NOT. found  .AND.  salsa )  THEN
-                CALL salsa_data_output_3d( av, do3d(av,ivar), found, local_pf )
+                CALL salsa_data_output_3d( av, do3d(av,ivar), found, local_pf, &
+                                           nzb_do, nzt_do )
                 resorted = .TRUE.
              ENDIF                  
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 3579)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 3582)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Move the control parameter "salsa" from salsa_mod.f90 to control_parameters
 ! 
 ! Former revisions:
@@ -165,5 +165,5 @@
                mask_j, mask_k, mask_size, mask_size_l, mask_start_l,           &
                mask_surface,                                                   &
-               max_masks, message_string, mid, nz_do3d, simulated_time
+               max_masks, message_string, mid, nz_do3d, salsa, simulated_time
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point
@@ -193,5 +193,5 @@
 
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_data_output_mask     
+        ONLY:  salsa_data_output_mask     
         
     USE surface_mod,                                                           &
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 3579)
+++ palm/trunk/SOURCE/header.f90	(revision 3582)
@@ -20,6 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
+!  
 ! Former revisions:
 ! -----------------
@@ -501,5 +502,5 @@
                
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_header        
+        ONLY:  salsa_header        
     
     USE spectra_mod,                                                           &
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 3579)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -655,5 +656,5 @@
              
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_init, salsa_init_arrays      
+        ONLY:  salsa_init, salsa_init_arrays      
     
     USE statistics,                                                            &
Index: palm/trunk/SOURCE/init_masks.f90
===================================================================
--- palm/trunk/SOURCE/init_masks.f90	(revision 3579)
+++ palm/trunk/SOURCE/init_masks.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -173,5 +174,5 @@
                mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
                mask_z_loop, max_masks,  message_string, mid,                   &
-               passive_scalar, ocean_mode, varnamelength
+               passive_scalar, ocean_mode, salsa, varnamelength
 
     USE grid_variables,                                                        &
@@ -190,10 +191,10 @@
 
     USE pegrid  
-
+    
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_check_data_output
         
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_check_data_output 
+        ONLY:  salsa_check_data_output 
 
     IMPLICIT NONE
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3579)
+++ palm/trunk/SOURCE/modules.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -1387,4 +1388,5 @@
     LOGICAL ::  run_control_header = .FALSE.                     !< onetime output of RUN_CONTROL header
     LOGICAL ::  run_coupled = .TRUE.                             !< internal switch telling PALM to run in coupled mode (i.e. to exchange surface data) in case of atmosphere-ocean coupling
+    LOGICAL ::  salsa = .FALSE.                                  !< switch for the sectional aerosol module salsa
     LOGICAL ::  scalar_rayleigh_damping = .TRUE.                 !< namelist parameter
     LOGICAL ::  sloping_surface = .FALSE.                        !< use sloped surface? (namelist parameter alpha_surface)
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3579)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -646,5 +647,5 @@
                message_string, mid, ntdim_2d_xy, ntdim_2d_xz,                  &
                ntdim_2d_yz, ntdim_3d, nz_do3d, ocean_mode, plant_canopy,       &
-               run_description_header, section, simulated_time,                &
+               run_description_header, salsa, section, simulated_time,         &
                simulated_time_at_begin, skip_time_data_output_av,              &
                skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
@@ -684,5 +685,5 @@
      
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_define_netcdf_grid             
+        ONLY:  salsa_define_netcdf_grid            
 
     USE spectra_mod,                                                           &
@@ -1147,5 +1148,5 @@
                       CALL salsa_define_netcdf_grid( domask(mid,av,i), found,  &
                                                      grid_x, grid_y, grid_z )
-                   ENDIF                  
+                   ENDIF                                   
 !
 !--                Now check for user-defined quantities
@@ -1910,5 +1911,5 @@
                                                      grid_y, grid_z )
                    ENDIF                 
-                   
+!                   
 !--                Check for user-defined quantities
                    IF ( .NOT. found )  THEN
@@ -2856,5 +2857,5 @@
                                                          grid_z )
                          ENDIF
-                         
+
 !
 !--                      Check for SALSA quantities
@@ -3758,5 +3759,4 @@
                                                             grid_z )
                       ENDIF
-                      
 !
 !--                   Check for SALSA quantities
@@ -4618,5 +4618,4 @@
                                                             grid_z )
                       ENDIF
-                      
 !
 !--                   Check for SALSA quantities
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3579)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3582)
@@ -22,5 +22,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Formatting adjustments
 ! 
 ! Former revisions:
@@ -1679,16 +1679,16 @@
                          ! to include also the dependecy to the radiation
                          ! in the plant canopy box
-                         pc_transpiration_rate(kk,j,i) =  - lsec                  &
-                                          * lad_s(kk,j,i) *                       &
-                                          SQRT( ( 0.5_wp * ( u(k,j,i) +           &
-                                                             u(k,j,i+1) )         &
-                                                )**2 +                            &
-                                                ( 0.5_wp * ( v(k,j,i) +           &
-                                                             v(k,j+1,i) )         &
-                                                )**2 +                            &
-                                                ( 0.5_wp * ( w(k-1,j,i) +         &
-                                                             w(k,j,i) )           &
-                                                )**2                              &
-                                              ) *                                 &
+                         pc_transpiration_rate(kk,j,i) =  - lsec               &
+                                          * lad_s(kk,j,i) *                    &
+                                          SQRT( ( 0.5_wp * ( u(k,j,i) +        &
+                                                             u(k,j,i+1) )      &
+                                                )**2 +                         &
+                                                ( 0.5_wp * ( v(k,j,i) +        &
+                                                             v(k,j+1,i) )      &
+                                                )**2 +                         &
+                                                ( 0.5_wp * ( w(k-1,j,i) +      &
+                                                             w(k,j,i) )        &
+                                                )**2                           &
+                                              ) *                              &
                                           ( q(k,j,i) - lsc )
                       ENDIF
@@ -2001,5 +2001,6 @@
                 DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
                    kk = k - k_wall  !- lad arrays are defined flat
-                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i)
+                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) -       &
+                                 pc_latent_rate(kk,j,i)
                 ENDDO
              ELSE
@@ -2024,16 +2025,16 @@
                    ! to include also the dependecy to the radiation
                    ! in the plant canopy box
-                   pc_transpiration_rate(kk,j,i) = - lsec                         &
-                                    * lad_s(kk,j,i) *                             &
-                                    SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
-                                                       u(k,j,i+1) )               &
-                                          )**2  +                                 &
-                                          ( 0.5_wp * ( v(k,j,i) +                 &
-                                                       v(k,j+1,i) )               &
-                                          )**2 +                                  &
-                                          ( 0.5_wp * ( w(k-1,j,i) +               &
-                                                       w(k,j,i) )                 &
-                                          )**2                                    &
-                                        ) *                                       &
+                   pc_transpiration_rate(kk,j,i) = - lsec                      &
+                                    * lad_s(kk,j,i) *                          &
+                                    SQRT( ( 0.5_wp * ( u(k,j,i) +              &
+                                                       u(k,j,i+1) )            &
+                                          )**2  +                              &
+                                          ( 0.5_wp * ( v(k,j,i) +              &
+                                                       v(k,j+1,i) )            &
+                                          )**2 +                               &
+                                          ( 0.5_wp * ( w(k-1,j,i) +            &
+                                                       w(k,j,i) )              &
+                                          )**2                                 &
+                                        ) *                                    &
                                     ( q(k,j,i) - lsc )
                 ENDIF
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 3579)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -373,5 +374,5 @@
                use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
                ws_scheme_sca, urban_surface, land_surface,                     &
-               time_since_reference_point
+               time_since_reference_point, salsa
 
     USE coriolis_mod,                                                          &
@@ -415,5 +416,5 @@
     USE salsa_mod,                                                             &
         ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time, nbins, &
-               ncc_tot, ngast, salsa, salsa_boundary_conds, salsa_diagnostics, &
+               ncc_tot, ngast, salsa_boundary_conds, salsa_diagnostics,        &
                salsa_driver, salsa_gas, salsa_gases_from_chem, salsa_tendency, &
                skip_time_do_salsa
@@ -548,5 +549,4 @@
 
     ENDIF       
-
 !
 !-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time 
Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3579)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -1087,5 +1088,5 @@
     
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_rrd_local
+        ONLY:  salsa_rrd_local
 
     USE surface_mod,                                                           &
@@ -1910,9 +1911,10 @@
                                            nyn_on_file, nysf, nysc,            &
                                            nys_on_file, found )
-!
-!--                Read salsa restart data
-                   IF ( .NOT. found  .AND.  salsa )  THEN
-                      CALL salsa_rrd_local
-                   ENDIF                                               
+
+                   IF ( .NOT. found  .AND.  salsa )  CALL salsa_rrd_local( i,  &
+                                           k, nxlf, nxlc, nxl_on_file, nxrf,   &
+                                           nxrc, nxr_on_file, nynf, nync,      & 
+                                           nyn_on_file, nysf, nysc,            &
+                                           nys_on_file, tmp_3d, found )                                        
 
 !
Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 3579)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 3582)
@@ -15,4 +15,5 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
+! Copyright 2018-2018 University of Helsinki
 ! Copyright 1997-2018 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
@@ -20,6 +21,10 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+! - Moved the control parameter "salsa" from salsa_mod.f90 to control_parameters
+! - Updated salsa_rrd_local and salsa_wrd_local
+! - Add target attribute
+! - Revise initialization in case of restarts
+! - Revise masked data output
+!
 ! Former revisions:
 ! -----------------
@@ -44,6 +49,6 @@
 ! Authors:
 ! --------
-! @author monakurppa
-!
+! @author Mona Kurppa (University of Helsinki)
+! 
 !
 ! Description:
@@ -60,6 +65,4 @@
 !>
 !> @todo Implement turbulent inflow of aerosols in inflow_turbulence.
-!> @todo Deposition on walls and horizontal surfaces calculated from the aerosol 
-!>       dry radius, not wet
 !> @todo Deposition on subgrid scale vegetation 
 !> @todo Deposition on vegetation calculated by default for deciduous broadleaf 
@@ -197,5 +200,4 @@
                                                 !< chemical components 
     LOGICAL ::  read_restart_data_salsa = .FALSE. !< read restart data for salsa
-    LOGICAL ::  salsa                 = .FALSE.   !< SALSA master switch
     LOGICAL ::  salsa_gases_from_chem = .FALSE.   !< Transfer the gaseous 
                                                   !< components to SALSA from  
@@ -420,9 +422,11 @@
                                                                     !< OC
 !-- Integrated:                                                                    
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  LDSA_av  !< lung deposited 
-                                                         !< surface area                                                    
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  Ntot_av  !< total number conc.
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  PM25_av  !< PM2.5
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  PM10_av  !< PM10
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  LDSA_av  !< lung-
+                                                                 !< deposited 
+                                                                 !< surface area                                                    
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  Ntot_av  !< total number
+                                                                 !< conc.
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  PM25_av  !< PM2.5
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  PM10_av  !< PM10
 !-- In the particle phase:    
     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_BC_av  !< black carbon 
@@ -435,9 +439,8 @@
     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_SS_av  !< sea salt
 !-- Bins:    
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  mbins_av  !< bin mass 
-                                                            !< concentration 
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  Nbins_av  !< bin number 
-                                                            !< concentration  
-        
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mbins_av !< bin mass  
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  Nbins_av !< bin number 
+
+    
 !
 !-- PALM interfaces:
@@ -564,4 +567,5 @@
            mconc_2, nbins, ncc, ncc_tot, nclim, nconc_2, ngast, prtcl, Ra_dry, &
            salsa_gas, sedim_vd
+           
 
  CONTAINS
@@ -688,6 +692,5 @@
 
 !
-!-- Set flag that indicates that the new module is switched on
-!-- Note that this parameter needs to be declared in modules.f90
+!-- Enable salsa (salsa switch in modules.f90)
     salsa = .TRUE.
 
@@ -1230,82 +1233,79 @@
     ENDDO
     
-    IF ( nldepo )  sedim_vd = 0.0_wp
+    IF ( nldepo )  sedim_vd = 0.0_wp  
+    
+    DO  b = 1, nbins
+       IF ( .NOT. read_restart_data_salsa )  aerosol_number(b)%conc = nclim
+       aerosol_number(b)%conc_p    = 0.0_wp
+       aerosol_number(b)%tconc_m   = 0.0_wp
+       aerosol_number(b)%flux_s    = 0.0_wp
+       aerosol_number(b)%diss_s    = 0.0_wp
+       aerosol_number(b)%flux_l    = 0.0_wp
+       aerosol_number(b)%diss_l    = 0.0_wp
+       aerosol_number(b)%init      = nclim
+       aerosol_number(b)%sums_ws_l = 0.0_wp
+    ENDDO
+    DO  c = 1, ncc_tot*nbins
+       IF ( .NOT. read_restart_data_salsa )  aerosol_mass(c)%conc = mclim
+       aerosol_mass(c)%conc_p    = 0.0_wp
+       aerosol_mass(c)%tconc_m   = 0.0_wp
+       aerosol_mass(c)%flux_s    = 0.0_wp
+       aerosol_mass(c)%diss_s    = 0.0_wp
+       aerosol_mass(c)%flux_l    = 0.0_wp
+       aerosol_mass(c)%diss_l    = 0.0_wp
+       aerosol_mass(c)%init      = mclim
+       aerosol_mass(c)%sums_ws_l = 0.0_wp
+    ENDDO
+    
+    IF ( .NOT. salsa_gases_from_chem )  THEN
+       DO  g = 1, ngast
+          salsa_gas(g)%conc_p    = 0.0_wp
+          salsa_gas(g)%tconc_m   = 0.0_wp
+          salsa_gas(g)%flux_s    = 0.0_wp
+          salsa_gas(g)%diss_s    = 0.0_wp
+          salsa_gas(g)%flux_l    = 0.0_wp
+          salsa_gas(g)%diss_l    = 0.0_wp
+          salsa_gas(g)%sums_ws_l = 0.0_wp
+       ENDDO
+       IF ( .NOT. read_restart_data_salsa )  THEN
+          salsa_gas(1)%conc = H2SO4_init
+          salsa_gas(2)%conc = HNO3_init
+          salsa_gas(3)%conc = NH3_init
+          salsa_gas(4)%conc = OCNV_init
+          salsa_gas(5)%conc = OCSV_init 
+       ENDIF
+!
+!--    Set initial value for gas compound tracers and initial values
+       salsa_gas(1)%init = H2SO4_init
+       salsa_gas(2)%init = HNO3_init
+       salsa_gas(3)%init = NH3_init
+       salsa_gas(4)%init = OCNV_init
+       salsa_gas(5)%init = OCSV_init      
+    ENDIF
+!
+!-- Aerosol radius in each bin: dry and wet (m)
+    Ra_dry = 1.0E-10_wp
 !    
-!-- Initilisation actions that are NOT conducted for restart runs
-    IF ( .NOT. read_restart_data_salsa )  THEN    
-    
-       DO  b = 1, nbins
-          aerosol_number(b)%conc      = nclim
-          aerosol_number(b)%conc_p    = 0.0_wp
-          aerosol_number(b)%tconc_m   = 0.0_wp
-          aerosol_number(b)%flux_s    = 0.0_wp
-          aerosol_number(b)%diss_s    = 0.0_wp
-          aerosol_number(b)%flux_l    = 0.0_wp
-          aerosol_number(b)%diss_l    = 0.0_wp
-          aerosol_number(b)%init      = nclim
-          aerosol_number(b)%sums_ws_l = 0.0_wp
+!-- Initialise aerosol tracers    
+    aero(:)%vhilim   = 0.0_wp
+    aero(:)%vlolim   = 0.0_wp
+    aero(:)%vratiohi = 0.0_wp
+    aero(:)%vratiolo = 0.0_wp
+    aero(:)%dmid     = 0.0_wp
+!
+!-- Initialise the sectional particle size distribution
+    CALL set_sizebins()
+!
+!-- Initialise location-dependent aerosol size distributions and 
+!-- chemical compositions:
+    CALL aerosol_init 
+!
+!-- Initalisation run of SALSA
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          CALL salsa_driver( i, j, 1 )
+          CALL salsa_diagnostics( i, j )
        ENDDO
-       DO  c = 1, ncc_tot*nbins
-          aerosol_mass(c)%conc      = mclim
-          aerosol_mass(c)%conc_p    = 0.0_wp
-          aerosol_mass(c)%tconc_m   = 0.0_wp
-          aerosol_mass(c)%flux_s    = 0.0_wp
-          aerosol_mass(c)%diss_s    = 0.0_wp
-          aerosol_mass(c)%flux_l    = 0.0_wp
-          aerosol_mass(c)%diss_l    = 0.0_wp
-          aerosol_mass(c)%init      = mclim
-          aerosol_mass(c)%sums_ws_l = 0.0_wp
-       ENDDO
-       
-       IF ( .NOT. salsa_gases_from_chem )  THEN
-          DO  g = 1, ngast
-             salsa_gas(g)%conc_p    = 0.0_wp
-             salsa_gas(g)%tconc_m   = 0.0_wp
-             salsa_gas(g)%flux_s    = 0.0_wp
-             salsa_gas(g)%diss_s    = 0.0_wp
-             salsa_gas(g)%flux_l    = 0.0_wp
-             salsa_gas(g)%diss_l    = 0.0_wp
-             salsa_gas(g)%sums_ws_l = 0.0_wp
-          ENDDO
-       
-!
-!--       Set initial value for gas compound tracers and initial values
-          salsa_gas(1)%conc = H2SO4_init
-          salsa_gas(1)%init = H2SO4_init
-          salsa_gas(2)%conc = HNO3_init
-          salsa_gas(2)%init = HNO3_init
-          salsa_gas(3)%conc = NH3_init
-          salsa_gas(3)%init = NH3_init
-          salsa_gas(4)%conc = OCNV_init
-          salsa_gas(4)%init = OCNV_init
-          salsa_gas(5)%conc = OCSV_init
-          salsa_gas(5)%init = OCSV_init      
-       ENDIF
-!
-!--    Aerosol radius in each bin: dry and wet (m)
-       Ra_dry = 1.0E-10_wp
-!    
-!--    Initialise aerosol tracers    
-       aero(:)%vhilim   = 0.0_wp
-       aero(:)%vlolim   = 0.0_wp
-       aero(:)%vratiohi = 0.0_wp
-       aero(:)%vratiolo = 0.0_wp
-       aero(:)%dmid     = 0.0_wp
-!
-!--    Initialise the sectional particle size distribution
-       CALL set_sizebins()
-!
-!--    Initialise location-dependent aerosol size distributions and 
-!--    chemical compositions:
-       CALL aerosol_init 
-!
-!--    Initalisation run of SALSA
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             CALL salsa_driver( i, j, 1 )
-             CALL salsa_diagnostics( i, j )
-          ENDDO
-       ENDDO 
-    ENDIF
+    ENDDO 
 !
 !-- Set the aerosol and gas sources
@@ -1419,6 +1419,6 @@
     
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  get_attribute, netcdf_data_input_get_dimension_length,          &
-               get_variable, open_read_file
+        ONLY:  get_attribute, get_variable,                                    &
+               netcdf_data_input_get_dimension_length, open_read_file
     
     IMPLICIT NONE
@@ -1485,5 +1485,6 @@
 !
 !--       Input heights   
-          CALL netcdf_data_input_get_dimension_length( id_faero, nz_file, "profile_z" ) 
+          CALL netcdf_data_input_get_dimension_length( id_faero, nz_file,      &
+                                                       "profile_z" ) 
           
           ALLOCATE( pr_z(nz_file), pr_mass_fracs_a(maxspec,nz_file),           &
@@ -1628,5 +1629,6 @@
 !
 !--       Input heights   
-          CALL netcdf_data_input_get_dimension_length( id_fchem, nz_file, "profile_z" ) 
+          CALL netcdf_data_input_get_dimension_length( id_fchem, nz_file,      &
+                                                       "profile_z" ) 
           ALLOCATE( pr_z(nz_file), pr_gas(ngast,nz_file) ) 
           CALL get_variable( id_fchem, 'profile_z', pr_z ) 
@@ -2025,5 +2027,7 @@
 !> This routine reads the respective restart data.
 !------------------------------------------------------------------------------!
- SUBROUTINE salsa_rrd_local  
+ SUBROUTINE salsa_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,        &
+                             nxr_on_file, nynf, nync, nyn_on_file, nysf,       &
+                             nysc, nys_on_file, tmp_3d, found )
 
     
@@ -2035,42 +2039,57 @@
     INTEGER(iwp) ::  g  !<
     INTEGER(iwp) ::  i  !<
-    INTEGER(iwp) ::  j  !<
-    INTEGER(iwp) ::  k  !<   
+    INTEGER(iwp) ::  k  !<
+    INTEGER(iwp) ::  nxlc            !<
+    INTEGER(iwp) ::  nxlf            !<
+    INTEGER(iwp) ::  nxl_on_file     !<
+    INTEGER(iwp) ::  nxrc            !<
+    INTEGER(iwp) ::  nxrf            !<
+    INTEGER(iwp) ::  nxr_on_file     !<
+    INTEGER(iwp) ::  nync            !<
+    INTEGER(iwp) ::  nynf            !<
+    INTEGER(iwp) ::  nyn_on_file     !<
+    INTEGER(iwp) ::  nysc            !<
+    INTEGER(iwp) ::  nysf            !<
+    INTEGER(iwp) ::  nys_on_file     !<
+
+    LOGICAL, INTENT(OUT)  ::  found
+
+    REAL(wp), &
+       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
+       
+    found = .FALSE.
     
     IF ( read_restart_data_salsa )  THEN
-       READ ( 13 )  field_char
-
-       DO  WHILE ( TRIM( field_char ) /= '*** end salsa ***' )
+    
+       SELECT CASE ( restart_string(1:length) )
        
-          DO b = 1, nbins
-             READ ( 13 )  aero(b)%vlolim
-             READ ( 13 )  aero(b)%vhilim
-             READ ( 13 )  aero(b)%dmid
-             READ ( 13 )  aero(b)%vratiohi
-             READ ( 13 )  aero(b)%vratiolo
-          ENDDO
-
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO k = nzb+1, nzt
-                   DO  b = 1, nbins
-                      READ ( 13 )  aerosol_number(b)%conc(k,j,i)
-                      DO  c = 1, ncc_tot
-                         READ ( 13 )  aerosol_mass((c-1)*nbins+b)%conc(k,j,i)
-                      ENDDO
-                   ENDDO
-                   IF ( .NOT. salsa_gases_from_chem )  THEN
-                      DO  g = 1, ngast
-                         READ ( 13 )  salsa_gas(g)%conc(k,j,i)
-                      ENDDO 
-                   ENDIF
-                ENDDO
+          CASE ( 'aerosol_number' )
+             DO  b = 1, nbins
+                IF ( k == 1 )  READ ( 13 ) tmp_3d
+                aerosol_number(b)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 
+                               tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                found = .TRUE.
              ENDDO
-          ENDDO
-
-          READ ( 13 )  field_char
-
-       ENDDO
        
+          CASE ( 'aerosol_mass' )
+             DO  c = 1, ncc_tot * nbins
+                IF ( k == 1 )  READ ( 13 ) tmp_3d
+                aerosol_mass(c)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 
+                               tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                found = .TRUE.
+             ENDDO
+          
+          CASE ( 'salsa_gas' )
+             DO  g = 1, ngast
+                IF ( k == 1 )  READ ( 13 ) tmp_3d
+                salsa_gas(g)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  & 
+                               tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                found = .TRUE.
+             ENDDO
+             
+          CASE DEFAULT
+             found = .FALSE.
+             
+       END SELECT
     ENDIF
 
@@ -2093,37 +2112,21 @@
     INTEGER(iwp) ::  c  !<
     INTEGER(iwp) ::  g  !<
-    INTEGER(iwp) ::  i  !<
-    INTEGER(iwp) ::  j  !<
-    INTEGER(iwp) ::  k  !<
     
     IF ( write_binary  .AND.  write_binary_salsa )  THEN
-       
-       DO b = 1, nbins
-          WRITE ( 14 )  aero(b)%vlolim
-          WRITE ( 14 )  aero(b)%vhilim
-          WRITE ( 14 )  aero(b)%dmid
-          WRITE ( 14 )  aero(b)%vratiohi
-          WRITE ( 14 )  aero(b)%vratiolo
+    
+       CALL wrd_write_string( 'aerosol_number' )
+       DO  b = 1, nbins
+          WRITE ( 14 )  aerosol_number(b)%conc
        ENDDO
        
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                DO  b = 1, nbins
-                   WRITE ( 14 )  aerosol_number(b)%conc(k,j,i)
-                   DO  c = 1, ncc_tot
-                      WRITE ( 14 )  aerosol_mass((c-1)*nbins+b)%conc(k,j,i)
-                   ENDDO
-                ENDDO
-                IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  g = 1, ngast
-                      WRITE ( 14 )  salsa_gas(g)%conc(k,j,i)
-                   ENDDO 
-                ENDIF
-             ENDDO
-          ENDDO
+       CALL wrd_write_string( 'aerosol_mass' )
+       DO  c = 1, nbins*ncc_tot
+          WRITE ( 14 )  aerosol_mass(c)%conc
        ENDDO
        
-       WRITE ( 14 )  '*** end salsa ***   '
+       CALL wrd_write_string( 'salsa_gas' )
+       DO  g = 1, ngast
+          WRITE ( 14 )  salsa_gas(g)%conc
+       ENDDO
           
     ENDIF
@@ -2765,33 +2768,33 @@
     USE arrays_3d,                                                             &
         ONLY: p, pt, q, zu
-        
+       
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  barometric_formula, exner_function, ideal_gas_law_rho, magnus 
-        
+        ONLY:  barometric_formula, exner_function, ideal_gas_law_rho, magnus
+       
     USE control_parameters,                                                    &
         ONLY: pt_surface, surface_pressure
-        
+       
     IMPLICIT NONE
    
     INTEGER(iwp), INTENT(in) ::  i
-    INTEGER(iwp), INTENT(in) ::  j 
+    INTEGER(iwp), INTENT(in) ::  j
     REAL(wp), DIMENSION(:), INTENT(inout) ::  adn_ij
-    REAL(wp), DIMENSION(:), INTENT(inout) ::  p_ij        
+    REAL(wp), DIMENSION(:), INTENT(inout) ::  p_ij       
     REAL(wp), DIMENSION(:), INTENT(inout) ::  temp_ij
     REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cw_ij
-    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cs_ij 
-    REAL(wp), DIMENSION(nzb:nzt+1) ::  e_s !< saturation vapour pressure 
+    REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::  cs_ij
+    REAL(wp), DIMENSION(nzb:nzt+1) ::  e_s !< saturation vapour pressure
                                            !< over water (Pa)
     REAL(wp) ::  t_surface !< absolute surface temperature (K)
 !
 !-- Pressure p_ijk (Pa) = hydrostatic pressure + perturbation pressure (p)
-    t_surface = pt_surface * exner_function( surface_pressure )
-    p_ij(:) = 100.0_wp * barometric_formula( zu, t_surface, surface_pressure ) &
+    t_surface = pt_surface * exner_function( surface_pressure * 100.0_wp )
+    p_ij(:) = barometric_formula( zu, t_surface, surface_pressure * 100.0_wp ) &
               + p(:,j,i)
 !             
 !-- Absolute ambient temperature (K)
-    temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) )        
-!
-!-- Air density 
+    temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) )       
+!
+!-- Air density
     adn_ij(:) = ideal_gas_law_rho( p_ij(:), temp_ij(:) )
 !
@@ -2803,8 +2806,9 @@
 !-- Saturation mixing ratio r_s (kg/kg) from vapour pressure at temp (Pa)
     IF ( PRESENT( cs_ij ) )  THEN
-       e_s(:) = magnus( temp_ij(:) ) 
-       cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:) - e_s(:) ) ) * adn_ij(:)  
+       e_s(:) = 611.0_wp * EXP( alv_d_rv * ( 3.6609E-3_wp - 1.0_wp /           &
+                temp_ij(:) ) )! magnus( temp_ij(:) )
+       cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:) - e_s(:) ) ) * adn_ij(:) 
     ENDIF
-    
+
  END SUBROUTINE salsa_thrm_ij 
 
@@ -7877,6 +7881,6 @@
     
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  get_attribute, netcdf_data_input_get_dimension_length,          &
-               get_variable, open_read_file
+        ONLY:  get_attribute, get_variable,                                    &
+               netcdf_data_input_get_dimension_length, open_read_file
     
     USE surface_mod,                                                           &
@@ -7985,5 +7989,6 @@
 !                              
 !--       Index of gaseous compounds
-          CALL netcdf_data_input_get_dimension_length( id_fchem, ng_file, "nspecies" )  
+          CALL netcdf_data_input_get_dimension_length( id_fchem, ng_file,      &
+                                                       "nspecies" )  
           IF ( ng_file < 5 )  THEN
              message_string = 'Some gaseous emissions missing.'
@@ -8200,5 +8205,6 @@
 !
 !--    Emission time step
-       CALL netcdf_data_input_get_dimension_length( id_faero, n_dt, 'dt_emission' ) 
+       CALL netcdf_data_input_get_dimension_length( id_faero, n_dt,            &
+                                                    'dt_emission' ) 
        IF ( n_dt > 1 )  THEN
           CALL location_message( '    salsa_set_source: hourly emission data'//&
@@ -8410,5 +8416,5 @@
        DO  n = 1, ncat_emission
           DO  g = 1, ngast
-             IF ( .NOT. salsa_gas(g)%source(n,j,i) > 0.0_wp )  THEN
+             IF ( salsa_gas(g)%source(n,j,i) < 0.0_wp )  THEN
                 salsa_gas(g)%source(n,j,i) = 0.0_wp
                 CYCLE
@@ -9138,6 +9144,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      PM25_av(k,j,i) = PM25_av(k,j,i)                          &
-                                        / REAL( average_count_3d, KIND=wp )
+                      PM25_av(k,j,i) = PM25_av(k,j,i) /                        &
+                                       REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -9148,6 +9154,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      PM10_av(k,j,i) = PM10_av(k,j,i)                          &
-                                        / REAL( average_count_3d, KIND=wp )
+                      PM10_av(k,j,i) = PM10_av(k,j,i) /                        &
+                                       REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -9168,6 +9174,6 @@
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i)         &
-                                             / REAL( average_count_3d, KIND=wp )
+                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /       &
+                                                 REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -9188,10 +9194,11 @@
 !> Subroutine defining 2D output variables
 !------------------------------------------------------------------------------!
- SUBROUTINE salsa_data_output_2d( av, variable, found, grid, mode,             &
-                                      local_pf, two_d )
+ SUBROUTINE salsa_data_output_2d( av, variable, found, grid, mode, local_pf,   &
+                                  two_d, nzb_do, nzt_do )
  
     USE indices
 
     USE kinds
+
 
     IMPLICIT NONE
@@ -9202,26 +9209,30 @@
     CHARACTER (LEN=5) ::  vari       !<  trimmed format of variable
 
-    INTEGER(iwp) ::  av   !< 
-    INTEGER(iwp) ::  b    !<
-    INTEGER(iwp) ::  c    !< 
-    INTEGER(iwp) ::  i    !<
-    INTEGER(iwp) ::  icc  !< index of a chemical compound
-    INTEGER(iwp) ::  j    !< 
-    INTEGER(iwp) ::  k    !< 
-
-    LOGICAL ::  found   !< 
-    LOGICAL ::  two_d   !< flag parameter that indicates 2D variables 
-                        !< (horizontal cross sections)
-    
-    REAL(wp) ::  df       !< For calculating LDSA: fraction of particles 
-                          !< depositing in the alveolar (or tracheobronchial)
-                          !< region of the lung. Depends on the particle size 
-    REAL(wp) ::  mean_d   !< Particle diameter in micrometres
-    REAL(wp) ::  nc       !< Particle number concentration in units 1/cm**3
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< local
-       !< array to which output data is resorted to
-    REAL(wp) ::  temp_bin !< 
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to 
-                                                     !< selected output variable
+    INTEGER(iwp) ::  av      !< 
+    INTEGER(iwp) ::  b       !< running index: size bins
+    INTEGER(iwp) ::  c       !< running index: mass bins
+    INTEGER(iwp) ::  i       !<
+    INTEGER(iwp) ::  icc     !< index of a chemical compound
+    INTEGER(iwp) ::  j       !< 
+    INTEGER(iwp) ::  k       !< 
+    INTEGER(iwp) ::  nzb_do  !< 
+    INTEGER(iwp) ::  nzt_do  !< 
+
+    LOGICAL ::  found        !< 
+    LOGICAL ::  two_d        !< flag parameter that indicates 2D variables 
+                             !< (horizontal cross sections)
+    
+    REAL(wp) ::  df          !< For calculating LDSA: fraction of particles 
+                             !< depositing in the alveolar (or tracheobronchial)
+                             !< region of the lung. Depends on the particle size 
+    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute                          
+    REAL(wp) ::  mean_d      !< Particle diameter in micrometres
+    REAL(wp) ::  nc          !< Particle number concentration in units 1/cm**3
+    REAL(wp) ::  temp_bin    !< temporary array for calculating output variables
+    
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< output
+    
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted           !< pointer
+    
     
     found = .TRUE.
@@ -9238,7 +9249,7 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
+                DO  k = nzb_do, nzt_do
                    local_pf(i,j,k) = MERGE( salsa_gas(icc)%conc(k,j,i),        &
-                                            REAL( -999.0_wp, KIND = wp ),      &
+                                            REAL( fill_value, KIND = wp ),     &
                                             BTEST( wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
@@ -9253,7 +9264,7 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
+                DO  k = nzb_do, nzt_do
                    local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),             &
-                                            REAL( -999.0_wp, KIND = wp ),      &
+                                            REAL( fill_value, KIND = wp ),     &
                                             BTEST( wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
@@ -9268,5 +9279,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
+                DO  k = nzb_do, nzt_do
                    temp_bin = 0.0_wp
                    DO  b = 1, nbins
@@ -9286,6 +9297,6 @@
                       temp_bin = temp_bin + pi * mean_d**2.0_wp * df * nc 
                    ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin,  REAL( -999.0_wp,        &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = & 
+                                            wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9294,7 +9305,8 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   local_pf(i,j,k) = MERGE( LDSA_av(k,j,i), REAL( -999.0_wp,   &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do
+                   local_pf(i,j,k) = MERGE( LDSA_av(k,j,i), REAL( fill_value,  &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9303,12 +9315,28 @@
 
        IF ( mode == 'xy' )  grid = 'zu'
-    
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin1' )  THEN
+       
+    ELSEIF ( TRIM( variable(1:5) ) == 'N_bin' )  THEN
+       
+       vari = TRIM( variable( 6:LEN( TRIM( variable ) ) - 3 ) )
+    
+       IF ( TRIM( vari ) == '1' ) b = 1
+       IF ( TRIM( vari ) == '2' ) b = 2
+       IF ( TRIM( vari ) == '3' ) b = 3
+       IF ( TRIM( vari ) == '4' ) b = 4
+       IF ( TRIM( vari ) == '5' ) b = 5
+       IF ( TRIM( vari ) == '6' ) b = 6
+       IF ( TRIM( vari ) == '7' ) b = 7
+       IF ( TRIM( vari ) == '8' ) b = 8
+       IF ( TRIM( vari ) == '9' ) b = 9
+       IF ( TRIM( vari ) == '10' ) b = 10
+       IF ( TRIM( vari ) == '11' ) b = 11
+       IF ( TRIM( vari ) == '12' ) b = 12
+       
        IF ( av == 0 )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(1)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
+                DO  k = nzb_do, nzt_do                      
+                   local_pf(i,j,k) = MERGE( aerosol_number(b)%conc(k,j,i),     &
+                                            REAL( fill_value, KIND = wp ),     &
                                             BTEST( wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
@@ -9318,7 +9346,7 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,1),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
+                DO  k = nzb_do, nzt_do                      
+                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,b),                 &
+                                            REAL( fill_value, KIND = wp ),     &
                                             BTEST( wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
@@ -9326,13 +9354,18 @@
           ENDDO
        ENDIF
-    
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin2' )  THEN
+       
+       IF ( mode == 'xy' )  grid = 'zu'
+    
+    ELSEIF ( TRIM( variable(1:4) ) == 'Ntot' )  THEN
        IF ( av == 0 )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(2)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do
+                   temp_bin = 0.0_wp
+                   DO  b = 1, nbins
+                      temp_bin = temp_bin + aerosol_number(b)%conc(k,j,i)
+                   ENDDO
+                   local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = &
+                                            wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9341,21 +9374,44 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,2),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do
+                   local_pf(i,j,k) = MERGE( Ntot_av(k,j,i), REAL( fill_value,  &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
           ENDDO
        ENDIF
+
+       IF ( mode == 'xy' )  grid = 'zu'
+    
+    
+    ELSEIF ( TRIM( variable(1:5) ) == 'm_bin' )  THEN
        
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin3' )  THEN
+       vari = TRIM( variable( 6:LEN( TRIM( variable ) ) - 3 ) )
+    
+       IF ( TRIM( vari ) == '1' ) b = 1
+       IF ( TRIM( vari ) == '2' ) b = 2
+       IF ( TRIM( vari ) == '3' ) b = 3
+       IF ( TRIM( vari ) == '4' ) b = 4
+       IF ( TRIM( vari ) == '5' ) b = 5
+       IF ( TRIM( vari ) == '6' ) b = 6
+       IF ( TRIM( vari ) == '7' ) b = 7
+       IF ( TRIM( vari ) == '8' ) b = 8
+       IF ( TRIM( vari ) == '9' ) b = 9
+       IF ( TRIM( vari ) == '10' ) b = 10
+       IF ( TRIM( vari ) == '11' ) b = 11
+       IF ( TRIM( vari ) == '12' ) b = 12
+       
        IF ( av == 0 )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(3)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do   
+                   temp_bin = 0.0_wp
+                   DO  c = b, ncc_tot * nbins, nbins
+                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
+                   ENDDO
+                   local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value,        &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) )
                 ENDDO
              ENDDO
@@ -9364,547 +9420,14 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,3),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do                      
+                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,b), REAL( fill_value,&
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
           ENDDO
        ENDIF
-    
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin4' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(4)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,4),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
        
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin5' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(5)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,5),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin6' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(6)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,6),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin7' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(7)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,7),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin8' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(8)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,8),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'N_bin9' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(9)%conc(k,j,i),     &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,9),                 &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-    
-    ELSEIF ( TRIM( variable(1:7) ) == 'N_bin10' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(10)%conc(k,j,i),    &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,10),                &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:7) ) == 'N_bin11' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(11)%conc(k,j,i),    &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,11),                &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:7) ) == 'N_bin12' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( aerosol_number(12)%conc(k,j,i),    &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,12),                &
-                                            REAL( -999.0_wp, KIND = wp ),      &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-    
-    ELSEIF ( TRIM( variable(1:4) ) == 'Ntot' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   temp_bin = 0.0_wp
-                   DO  b = 1, nbins
-                      temp_bin = temp_bin + aerosol_number(b)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   local_pf(i,j,k) = MERGE( Ntot_av(k,j,i), REAL( -999.0_wp,   &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-
        IF ( mode == 'xy' )  grid = 'zu'
-    
-    
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin1' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 1, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,1), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-    
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin2' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 2, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,2), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin3' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 3, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,3), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin4' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 4, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,4), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin5' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 5, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,5), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin6' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 6, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,6), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin7' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 7, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,7), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin8' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 8, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,8), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:6) ) == 'm_bin9' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 9, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,9), REAL( -999.0_wp,&
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:7) ) == 'm_bin10' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 10, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,10), REAL(          &
-                       -999.0_wp, KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:7) ) == 'm_bin11' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 11, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,11), REAL(          &
-                       -999.0_wp, KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-       
-    ELSEIF ( TRIM( variable(1:7) ) == 'm_bin12' )  THEN
-       IF ( av == 0 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1   
-                   temp_bin = 0.0_wp
-                   DO  c = 12, ncc_tot*nbins, nbins
-                      temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                   ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSE
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1                      
-                   local_pf(i,j,k) = MERGE( mbins_av(k,j,i,12), REAL(          &
-                       -999.0_wp, KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
     
     ELSEIF ( TRIM( variable(1:5) ) == 'PM2.5' )  THEN
@@ -9912,5 +9435,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
+                DO  k = nzb_do, nzt_do
                    temp_bin = 0.0_wp
                    DO  b = 1, nbins
@@ -9921,6 +9444,7 @@
                       ENDIF
                    ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,         &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value,        &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9929,7 +9453,8 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   local_pf(i,j,k) = MERGE( PM25_av(k,j,i), REAL( -999.0_wp,   &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do
+                   local_pf(i,j,k) = MERGE( PM25_av(k,j,i), REAL( fill_value,  &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9944,5 +9469,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
+                DO  k = nzb_do, nzt_do
                    temp_bin = 0.0_wp
                    DO  b = 1, nbins
@@ -9953,6 +9478,7 @@
                       ENDIF
                    ENDDO
-                   local_pf(i,j,k) = MERGE( temp_bin,  REAL( -999.0_wp,        &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   local_pf(i,j,k) = MERGE( temp_bin,  REAL( fill_value,       &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9961,7 +9487,8 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   local_pf(i,j,k) = MERGE( PM10_av(k,j,i), REAL( -999.0_wp,   &
-                                 KIND = wp ),  BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                DO  k = nzb_do, nzt_do
+                   local_pf(i,j,k) = MERGE( PM10_av(k,j,i), REAL( fill_value,  &
+                                            KIND = wp ), BTEST(                &
+                                            wall_flags_0(k,j,i), 0 ) ) 
                 ENDDO
              ENDDO
@@ -9978,11 +9505,12 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       temp_bin = 0.0_wp
                       DO  c = ( icc-1 )*nbins+1, icc*nbins, 1
                          temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
                       ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin, REAL( -999.0_wp,      &
-                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                      local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value,     &
+                                               KIND = wp ), BTEST(             &
+                                               wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
                 ENDDO
@@ -9998,7 +9526,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),          &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10007,5 +9535,5 @@
           ENDIF
        ELSE
-          local_pf = 0.0_wp 
+          local_pf = fill_value
        ENDIF
 
@@ -10027,9 +9555,11 @@
 !> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
- SUBROUTINE salsa_data_output_3d( av, variable, found, local_pf )
+ SUBROUTINE salsa_data_output_3d( av, variable, found, local_pf, nzb_do,       &
+                                  nzt_do )
 
     USE indices
 
     USE kinds
+    
 
     IMPLICIT NONE
@@ -10037,24 +9567,28 @@
     CHARACTER (LEN=*), INTENT(in) ::  variable   !< 
     
-    INTEGER(iwp) ::  av   !< 
-    INTEGER(iwp) ::  c    !<
-    INTEGER(iwp) ::  i    !<
-    INTEGER(iwp) ::  icc  !< index of a chemical compound
-    INTEGER(iwp) ::  j    !< 
-    INTEGER(iwp) ::  k    !< 
-    INTEGER(iwp) ::  n    !<
-
-    LOGICAL ::  found   !< 
-    REAL(wp) ::  df       !< For calculating LDSA: fraction of particles 
-                          !< depositing in the alveolar (or tracheobronchial)
-                          !< region of the lung. Depends on the particle size 
-    REAL(wp) ::  mean_d   !< Particle diameter in micrometres
-    REAL(wp) ::  nc       !< Particle number concentration in units 1/cm**3
-
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf  !< local
-                                  !< array to which output data is resorted to
-    REAL(wp) ::  temp_bin  !< 
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to 
-                                                     !< selected output variable
+    INTEGER(iwp) ::  av      !<
+    INTEGER(iwp) ::  b       !< running index: size bins    
+    INTEGER(iwp) ::  c       !< running index: mass bins
+    INTEGER(iwp) ::  i       !<
+    INTEGER(iwp) ::  icc     !< index of a chemical compound
+    INTEGER(iwp) ::  j       !< 
+    INTEGER(iwp) ::  k       !< 
+    INTEGER(iwp) ::  nzb_do  !< 
+    INTEGER(iwp) ::  nzt_do  !<   
+
+    LOGICAL ::  found      !< 
+    
+    REAL(wp) ::  df        !< For calculating LDSA: fraction of particles 
+                           !< depositing in the alveolar (or tracheobronchial)
+                           !< region of the lung. Depends on the particle size 
+    REAL(wp) ::  fill_value = -9999.0_wp   !< value for the _FillValue attribute
+    REAL(wp) ::  mean_d    !< Particle diameter in micrometres
+    REAL(wp) ::  nc        !< Particle number concentration in units 1/cm**3
+    REAL(wp) ::  temp_bin  !< temporary array for calculating output variables    
+
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< local 
+    
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< pointer
+                                                     
        
     found     = .TRUE.
@@ -10073,7 +9607,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( salsa_gas(icc)%conc(k,j,i),     &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10088,7 +9622,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),          &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10101,19 +9635,19 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
+                      DO  b = 1, nbins
 !                      
 !--                      Diameter in micrometres
-                         mean_d = 1.0E+6_wp * Ra_dry(k,j,i,n) * 2.0_wp 
+                         mean_d = 1.0E+6_wp * Ra_dry(k,j,i,b) * 2.0_wp 
 !                               
 !--                      Deposition factor: alveolar                             
                          df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp *     &
-                                ( LOG( mean_d ) + 2.84_wp )**2.0_wp )          &
-                                  + 19.11_wp * EXP( -0.482_wp *                &
-                                  ( LOG( mean_d ) - 1.362_wp )**2.0_wp ) )
+                              ( LOG( mean_d ) + 2.84_wp )**2.0_wp ) + 19.11_wp &
+                              * EXP( -0.482_wp * ( LOG( mean_d ) - 1.362_wp    &
+                                )**2.0_wp ) )
 !                                    
 !--                      Number concentration in 1/cm3
-                         nc = 1.0E-6_wp * aerosol_number(n)%conc(k,j,i)
+                         nc = 1.0E-6_wp * aerosol_number(b)%conc(k,j,i)
 !                         
 !--                      Lung-deposited surface area LDSA (units mum2/cm3)
@@ -10121,5 +9655,5 @@
                       ENDDO
                       local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10129,7 +9663,44 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( LDSA_av(k,j,i),                 &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          
+       CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4',   'N_bin5',  'N_bin6',   &
+              'N_bin7', 'N_bin8', 'N_bin9', 'N_bin10' , 'N_bin11', 'N_bin12' )
+          IF ( TRIM( variable(6:) ) == '1' ) b = 1
+          IF ( TRIM( variable(6:) ) == '2' ) b = 2
+          IF ( TRIM( variable(6:) ) == '3' ) b = 3
+          IF ( TRIM( variable(6:) ) == '4' ) b = 4
+          IF ( TRIM( variable(6:) ) == '5' ) b = 5
+          IF ( TRIM( variable(6:) ) == '6' ) b = 6
+          IF ( TRIM( variable(6:) ) == '7' ) b = 7
+          IF ( TRIM( variable(6:) ) == '8' ) b = 8
+          IF ( TRIM( variable(6:) ) == '9' ) b = 9
+          IF ( TRIM( variable(6:) ) == '10' ) b = 10
+          IF ( TRIM( variable(6:) ) == '11' ) b = 11
+          IF ( TRIM( variable(6:) ) == '12' ) b = 12
+          
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do                      
+                      local_pf(i,j,k) = MERGE( aerosol_number(b)%conc(k,j,i),  &
+                                               REAL( fill_value, KIND = wp ),  &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do                      
+                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,b),              &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10142,11 +9713,11 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins                         
-                         temp_bin = temp_bin + aerosol_number(n)%conc(k,j,i)
+                      DO  b = 1, nbins                         
+                         temp_bin = temp_bin + aerosol_number(b)%conc(k,j,i)
                       ENDDO
                       local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10156,7 +9727,48 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( Ntot_av(k,j,i),                 &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          
+       CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4',   'm_bin5',  'm_bin6',   &
+              'm_bin7', 'm_bin8', 'm_bin9', 'm_bin10' , 'm_bin11', 'm_bin12' )
+          IF ( TRIM( variable(6:) ) == '1' ) b = 1
+          IF ( TRIM( variable(6:) ) == '2' ) b = 2
+          IF ( TRIM( variable(6:) ) == '3' ) b = 3
+          IF ( TRIM( variable(6:) ) == '4' ) b = 4
+          IF ( TRIM( variable(6:) ) == '5' ) b = 5
+          IF ( TRIM( variable(6:) ) == '6' ) b = 6
+          IF ( TRIM( variable(6:) ) == '7' ) b = 7
+          IF ( TRIM( variable(6:) ) == '8' ) b = 8
+          IF ( TRIM( variable(6:) ) == '9' ) b = 9
+          IF ( TRIM( variable(6:) ) == '10' ) b = 10
+          IF ( TRIM( variable(6:) ) == '11' ) b = 11
+          IF ( TRIM( variable(6:) ) == '12' ) b = 12
+          
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do   
+                      temp_bin = 0.0_wp
+                      DO  c = b, ncc_tot * nbins, nbins
+                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
+                      ENDDO
+                      local_pf(i,j,k) = MERGE( temp_bin,                       &
+                                               REAL( fill_value, KIND = wp ),  &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do                      
+                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,b),              &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10169,9 +9781,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
-                         IF ( 2.0_wp * Ra_dry(k,j,i,n) <= 2.5E-6_wp )  THEN
-                            DO  c = n, nbins*ncc, nbins
+                      DO  b = 1, nbins
+                         IF ( 2.0_wp * Ra_dry(k,j,i,b) <= 2.5E-6_wp )  THEN
+                            DO  c = b, nbins * ncc, nbins
                                temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
                             ENDDO
@@ -10179,5 +9791,5 @@
                       ENDDO
                       local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10187,7 +9799,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( PM25_av(k,j,i),                 &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10200,9 +9812,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
-                         IF ( 2.0_wp * Ra_dry(k,j,i,n) <= 10.0E-6_wp )  THEN
-                            DO  c = n, nbins*ncc, nbins
+                      DO  b = 1, nbins
+                         IF ( 2.0_wp * Ra_dry(k,j,i,b) <= 10.0E-6_wp )  THEN
+                            DO  c = b, nbins * ncc, nbins
                                temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
                             ENDDO
@@ -10210,5 +9822,5 @@
                       ENDDO
                       local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10218,607 +9830,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = MERGE( PM10_av(k,j,i),                 &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin1' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(1)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,1),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin2' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1 
-                      local_pf(i,j,k) = MERGE( aerosol_number(2)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,2),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin3' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(3)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,3),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin4' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1    
-                      local_pf(i,j,k) = MERGE( aerosol_number(4)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,4),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin5' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(5)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,5),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin6' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(6)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,6),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin7' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(7)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,7),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin8' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                  
-                      local_pf(i,j,k) = MERGE( aerosol_number(8)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,8),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin9' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(9)%conc(k,j,i),  &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,9),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin10' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(10)%conc(k,j,i), &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,10),             &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin11' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(11)%conc(k,j,i), &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,11),             &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin12' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( aerosol_number(12)%conc(k,j,i), &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( Nbins_av(k,j,i,12),             &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin1' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 1, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,1),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin2' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 2, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,2),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin3' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 3, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,3),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin4' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 4, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,4),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin5' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 5, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,5),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin6' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 6, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,6),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin7' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 7, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,7),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin8' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 8, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,8),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin9' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 9, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,9),              &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin10' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 10, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,10),             &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin11' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 11, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,11),             &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) ) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin12' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1   
-                      temp_bin = 0.0_wp
-                      DO  c = 12, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
-                      ENDDO
-                      local_pf(i,j,k) = MERGE( temp_bin,                       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb, nzt+1                      
-                      local_pf(i,j,k) = MERGE( mbins_av(k,j,i,12),             &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                    ENDDO
@@ -10833,5 +9845,5 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          temp_bin = 0.0_wp
                          DO  c = ( icc-1 )*nbins+1, icc*nbins                          
@@ -10839,5 +9851,5 @@
                          ENDDO
                          local_pf(i,j,k) = MERGE( temp_bin,                    &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                       ENDDO
@@ -10854,7 +9866,7 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      DO  k = nzb, nzt+1                      
+                      DO  k = nzb_do, nzt_do                      
                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),       &
-                                               REAL( -999.0_wp, KIND = wp ),   &
+                                               REAL( fill_value, KIND = wp ),  &
                                                BTEST( wall_flags_0(k,j,i), 0 ) ) 
                       ENDDO
@@ -10878,34 +9890,45 @@
  SUBROUTINE salsa_data_output_mask( av, variable, found, local_pf )
  
+    USE arrays_3d,                                                             &
+        ONLY:  tend
+ 
     USE control_parameters,                                                    &
-        ONLY:  mask_size_l, mid
+        ONLY:  mask_size_l, mask_surface, mid
+        
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index_ji       
  
     IMPLICIT NONE
     
-    CHARACTER (LEN=*) ::  variable   !< 
-
-    INTEGER(iwp) ::  av   !< 
-    INTEGER(iwp) ::  c    !<
-    INTEGER(iwp) ::  i    !<
-    INTEGER(iwp) ::  icc  !< index of a chemical compound
-    INTEGER(iwp) ::  j    !< 
-    INTEGER(iwp) ::  k    !< 
-    INTEGER(iwp) ::  n    !<
-
-    LOGICAL  ::  found    !< 
-    REAL(wp) ::  df       !< For calculating LDSA: fraction of particles 
-                          !< depositing in the alveolar (or tracheobronchial)
-                          !< region of the lung. Depends on the particle size 
-    REAL(wp) ::  mean_d       !< Particle diameter in micrometres
-    REAL(wp) ::  nc       !< Particle number concentration in units 1/cm**3
-
-    REAL(wp),                                                                  &
-       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
-          local_pf   !< 
-    REAL(wp) ::  temp_bin   !< 
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to 
-                                                     !< selected output variable
+    CHARACTER(LEN=5) ::  grid      !< flag to distinquish between staggered grid
+    CHARACTER(LEN=*) ::  variable  !< 
+    CHARACTER(LEN=7) ::  vari      !< trimmed format of variable
+
+    INTEGER(iwp) ::  av              !< 
+    INTEGER(iwp) ::  b               !< loop index for aerosol size number bins
+    INTEGER(iwp) ::  c               !< loop index for chemical components
+    INTEGER(iwp) ::  i               !< loop index in x-direction
+    INTEGER(iwp) ::  icc             !< index of a chemical compound
+    INTEGER(iwp) ::  j               !< loop index in y-direction
+    INTEGER(iwp) ::  k               !< loop index in z-direction
+    INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
+    
+    LOGICAL ::  found      !< 
+    LOGICAL ::  resorted   !<
+    
+    REAL(wp) ::  df        !< For calculating LDSA: fraction of particles 
+                           !< depositing in the alveolar (or tracheobronchial)
+                           !< region of the lung. Depends on the particle size 
+    REAL(wp) ::  mean_d    !< Particle diameter in micrometres
+    REAL(wp) ::  nc        !< Particle number concentration in units 1/cm**3
+    REAL(wp) ::  temp_bin  !< temporary array for calculating output variables    
+
+    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !< 
+    
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< pointer
 
     found     = .TRUE.
+    resorted  = .FALSE.
+    grid      = 's'
     temp_bin  = 0.0_wp
 
@@ -10913,745 +9936,356 @@
     
        CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
+          vari = TRIM( variable )
           IF ( av == 0 )  THEN
-             IF ( TRIM( variable ) == 'g_H2SO4')  icc = 1
-             IF ( TRIM( variable ) == 'g_HNO3')   icc = 2
-             IF ( TRIM( variable ) == 'g_NH3')    icc = 3
-             IF ( TRIM( variable ) == 'g_OCNV')   icc = 4
-             IF ( TRIM( variable ) == 'g_OCSV')   icc = 5
-             
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      local_pf(i,j,k) = salsa_gas(icc)%conc(mask_k(mid,k),     &
-                                                    mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
+             IF ( vari == 'g_H2SO4')  to_be_resorted => salsa_gas(1)%conc
+             IF ( vari == 'g_HNO3')   to_be_resorted => salsa_gas(2)%conc
+             IF ( vari == 'g_NH3')    to_be_resorted => salsa_gas(3)%conc
+             IF ( vari == 'g_OCNV')   to_be_resorted => salsa_gas(4)%conc
+             IF ( vari == 'g_OCSV')   to_be_resorted => salsa_gas(5)%conc 
           ELSE
-             IF ( TRIM( variable(3:) ) == 'H2SO4' ) to_be_resorted => g_H2SO4_av
-             IF ( TRIM( variable(3:) ) == 'HNO3' )  to_be_resorted => g_HNO3_av   
-             IF ( TRIM( variable(3:) ) == 'NH3' )   to_be_resorted => g_NH3_av   
-             IF ( TRIM( variable(3:) ) == 'OCNV' )  to_be_resorted => g_OCNV_av   
-             IF ( TRIM( variable(3:) ) == 'OCSV' )  to_be_resorted => g_OCSV_av 
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      local_pf(i,j,k) = to_be_resorted(mask_k(mid,k),          &
-                                                    mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
+             IF ( vari == 'g_H2SO4') to_be_resorted => g_H2SO4_av
+             IF ( vari == 'g_HNO3')  to_be_resorted => g_HNO3_av   
+             IF ( vari == 'g_NH3')   to_be_resorted => g_NH3_av   
+             IF ( vari == 'g_OCNV')  to_be_resorted => g_OCNV_av   
+             IF ( vari == 'g_OCSV')  to_be_resorted => g_OCSV_av
           ENDIF
-       
+          
        CASE ( 'LDSA' )
           IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nz_do3d
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
+                      DO  b = 1, nbins
 !                      
 !--                      Diameter in micrometres
-                         mean_d = 1.0E+6_wp * Ra_dry(mask_k(mid,k),            &
-                                       mask_j(mid,j),mask_i(mid,i),n) * 2.0_wp
+                         mean_d = 1.0E+6_wp * Ra_dry(k,j,i,b) * 2.0_wp
 !                               
-!--                      Deposition factor: alveolar (use Ra_dry for the size??)                                
+!--                      Deposition factor: alveolar                                
                          df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp *     &
-                                ( LOG( mean_d ) + 2.84_wp )**2.0_wp )          &
-                                  + 19.11_wp * EXP( -0.482_wp *                &
-                                  ( LOG( mean_d ) - 1.362_wp )**2.0_wp ) )
+                              ( LOG( mean_d ) + 2.84_wp )**2.0_wp ) + 19.11_wp &
+                              * EXP( -0.482_wp * ( LOG( mean_d ) - 1.362_wp    &
+                                )**2.0_wp ) )
 !                                    
 !--                      Number concentration in 1/cm3
-                         nc = 1.0E-6_wp * aerosol_number(n)%conc(mask_k(mid,k),&
-                                                    mask_j(mid,j),mask_i(mid,i))
+                         nc = 1.0E-6_wp * aerosol_number(b)%conc(k,j,i)
 !                         
 !--                      Lung-deposited surface area LDSA (units mum2/cm3)
                          temp_bin = temp_bin + pi * mean_d**2.0_wp * df * nc 
                       ENDDO
-                      local_pf(i,j,k) = temp_bin
+                      tend(k,j,i) = temp_bin
                    ENDDO
                 ENDDO
              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = LDSA_av(mask_k(mid,k),                &
-                                                 mask_j(mid,j),mask_i(mid,i))
+             IF ( .NOT. mask_surface(mid) )  THEN   
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k),  mask_j(mid,j),&
+                                                 mask_i(mid,i) )
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
+             ELSE 
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k),&
+                                                      nzt+1 ),                 &
+                                                 mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+          ELSE
+             to_be_resorted => LDSA_av
           ENDIF
-       
+          
+       CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4',   'N_bin5',  'N_bin6',   &
+              'N_bin7', 'N_bin8', 'N_bin9', 'N_bin10' , 'N_bin11', 'N_bin12' )
+          IF ( TRIM( variable(6:) ) == '1' ) b = 1
+          IF ( TRIM( variable(6:) ) == '2' ) b = 2
+          IF ( TRIM( variable(6:) ) == '3' ) b = 3
+          IF ( TRIM( variable(6:) ) == '4' ) b = 4
+          IF ( TRIM( variable(6:) ) == '5' ) b = 5
+          IF ( TRIM( variable(6:) ) == '6' ) b = 6
+          IF ( TRIM( variable(6:) ) == '7' ) b = 7
+          IF ( TRIM( variable(6:) ) == '8' ) b = 8
+          IF ( TRIM( variable(6:) ) == '9' ) b = 9
+          IF ( TRIM( variable(6:) ) == '10' ) b = 10
+          IF ( TRIM( variable(6:) ) == '11' ) b = 11
+          IF ( TRIM( variable(6:) ) == '12' ) b = 12
+          
+          IF ( av == 0 )  THEN
+             IF ( .NOT. mask_surface(mid) )  THEN   
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = aerosol_number(b)%conc( mask_k(mid,k),&
+                                                                   mask_j(mid,j),&
+                                                                   mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE 
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = aerosol_number(b)%conc(             &
+                                           MIN( topo_top_ind+mask_k(mid,k),    &
+                                                nzt+1 ),                 &
+                                           mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+          ELSE
+             to_be_resorted => Nbins_av(:,:,:,b)
+          ENDIF
+          
        CASE ( 'Ntot' )
           IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nz_do3d
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
-                         temp_bin = temp_bin + aerosol_number(n)%conc(         &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                      DO  b = 1, nbins
+                         temp_bin = temp_bin + aerosol_number(b)%conc(k,j,i)
                       ENDDO
-                      local_pf(i,j,k) = temp_bin
+                      tend(k,j,i) = temp_bin
                    ENDDO
                 ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Ntot_av(mask_k(mid,k),                &
-                                                 mask_j(mid,j),mask_i(mid,i))
+             ENDDO 
+             IF ( .NOT. mask_surface(mid) )  THEN   
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k),  mask_j(mid,j),&
+                                                 mask_i(mid,i) )
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
+             ELSE 
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k),&
+                                                      nzt+1 ),                 &
+                                                 mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+          ELSE
+             to_be_resorted => Ntot_av
+          ENDIF
+          
+       CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4',   'm_bin5',  'm_bin6',   &
+              'm_bin7', 'm_bin8', 'm_bin9', 'm_bin10' , 'm_bin11', 'm_bin12' )
+          IF ( TRIM( variable(6:) ) == '1' ) b = 1
+          IF ( TRIM( variable(6:) ) == '2' ) b = 2
+          IF ( TRIM( variable(6:) ) == '3' ) b = 3
+          IF ( TRIM( variable(6:) ) == '4' ) b = 4
+          IF ( TRIM( variable(6:) ) == '5' ) b = 5
+          IF ( TRIM( variable(6:) ) == '6' ) b = 6
+          IF ( TRIM( variable(6:) ) == '7' ) b = 7
+          IF ( TRIM( variable(6:) ) == '8' ) b = 8
+          IF ( TRIM( variable(6:) ) == '9' ) b = 9
+          IF ( TRIM( variable(6:) ) == '10' ) b = 10
+          IF ( TRIM( variable(6:) ) == '11' ) b = 11
+          IF ( TRIM( variable(6:) ) == '12' ) b = 12
+          
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nz_do3d
+                      temp_bin = 0.0_wp
+                      DO  c = b, ncc_tot*nbins, nbins
+                         temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
+                      ENDDO
+                      tend(k,j,i) = temp_bin
+                   ENDDO
+                ENDDO
+             ENDDO   
+             IF ( .NOT. mask_surface(mid) )  THEN   
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k),  mask_j(mid,j),&
+                                                 mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE 
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k),&
+                                                      nzt+1 ),                 &
+                                                 mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+          ELSE
+             to_be_resorted => mbins_av(:,:,:,b)
           ENDIF
        
        CASE ( 'PM2.5' )
           IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nz_do3d
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
-                         IF ( 2.0_wp * Ra_dry(mask_k(mid,k),mask_j(mid,j),     &
-                              mask_i(mid,i),n) <= 2.5E-6_wp )  THEN
-                            DO  c = n, nbins*ncc, nbins
-                               temp_bin = temp_bin + aerosol_mass(c)%conc(     &
-                                     mask_k(mid,k), mask_j(mid,j),mask_i(mid,i))
+                      DO  b = 1, nbins
+                         IF ( 2.0_wp * Ra_dry(k,j,i,b) <= 2.5E-6_wp )  THEN
+                            DO  c = b, nbins * ncc, nbins
+                               temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
                             ENDDO
                          ENDIF
                       ENDDO
-                      local_pf(i,j,k) = temp_bin
+                      tend(k,j,i) = temp_bin
                    ENDDO
                 ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = PM25_av(mask_k(mid,k),                &
-                                                 mask_j(mid,j),mask_i(mid,i))
+             ENDDO 
+             IF ( .NOT. mask_surface(mid) )  THEN   
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k),  mask_j(mid,j),&
+                                                 mask_i(mid,i) )
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
+             ELSE 
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k),&
+                                                      nzt+1 ),                 &
+                                                 mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+          ELSE
+             to_be_resorted => PM25_av
           ENDIF
-       
+          
        CASE ( 'PM10' )
           IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nz_do3d
                       temp_bin = 0.0_wp
-                      DO  n = 1, nbins
-                         IF ( 2.0_wp * Ra_dry(mask_k(mid,k),mask_j(mid,j),     &
-                              mask_i(mid,i),n) <= 10.0E-6_wp )  THEN
-                            DO  c = n, nbins*ncc, nbins
-                               temp_bin = temp_bin + aerosol_mass(c)%conc(     &
-                                      mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                      DO  b = 1, nbins
+                         IF ( 2.0_wp * Ra_dry(k,j,i,b) <= 10.0E-6_wp )  THEN
+                            DO  c = b, nbins * ncc, nbins
+                               temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
                             ENDDO
                          ENDIF
                       ENDDO
-                      local_pf(i,j,k) = temp_bin
+                      tend(k,j,i) = temp_bin
                    ENDDO
                 ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = PM10_av(mask_k(mid,k),                &
-                                                 mask_j(mid,j),mask_i(mid,i))
+             ENDDO 
+             IF ( .NOT. mask_surface(mid) )  THEN   
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k),  mask_j(mid,j),&
+                                                 mask_i(mid,i) )
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin1' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(1)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i))
+             ELSE 
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k),&
+                                                      nzt+1 ),                 &
+                                                 mask_j(mid,j), mask_i(mid,i) )
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
+             ENDIF
+             resorted = .TRUE.
           ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),1)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin2' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(2)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),2)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin3' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(3)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),3)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin4' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(4)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),4)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin5' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(5)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),5)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin6' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(6)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),6)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin7' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(7)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),7)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin8' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(8)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),8)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin9' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(9)%conc(mask_k(mid,k),  &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),9)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin10' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(10)%conc(mask_k(mid,k), &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                 mask_j(mid,j),mask_i(mid,i),10)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'N_bin11' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(11)%conc(mask_k(mid,k), &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                 mask_j(mid,j),mask_i(mid,i),11)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'N_bin12' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)                      
-                      local_pf(i,j,k) = aerosol_number(12)%conc(mask_k(mid,k), &
-                                                 mask_j(mid,j),mask_i(mid,i)) 
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = Nbins_av(mask_k(mid,k),               &
-                                                 mask_j(mid,j),mask_i(mid,i),12)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin1' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 1, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),1)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin2' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 2, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),2)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin3' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 3, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),3)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin4' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 4, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),4)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin5' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 5, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),5)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin6' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 6, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),6)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin7' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 7, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),7)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin8' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 8, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),8)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin9' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 9, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                  mask_j(mid,j),mask_i(mid,i),9)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       
-       CASE ( 'm_bin10' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 10, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                 mask_j(mid,j),mask_i(mid,i),10)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin11' )
-         IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 11, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                 mask_j(mid,j),mask_i(mid,i),11)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-          
-       CASE ( 'm_bin12' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                      temp_bin = 0.0_wp
-                      DO  c = 12, ncc_tot*nbins, nbins
-                         temp_bin = temp_bin + aerosol_mass(c)%conc(           &
-                                    mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                      local_pf(i,j,k) = temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = mbins_av(mask_k(mid,k),               &
-                                                 mask_j(mid,j),mask_i(mid,i),12)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => PM10_av
           ENDIF
           
        CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
-          IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
-             icc = get_index( prtcl, TRIM( variable(3:) ) )
-             IF ( av == 0 )  THEN
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2)
-                      DO  k = 1, mask_size_l(mid,3)
+          IF ( av == 0 )  THEN
+             IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
+                icc = get_index( prtcl, TRIM( variable(3:) ) )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nz_do3d
                          temp_bin = 0.0_wp
                          DO  c = ( icc-1 )*nbins+1, icc*nbins 
-                            temp_bin = temp_bin + aerosol_mass(c)%conc(        &
-                                      mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                            temp_bin = temp_bin + aerosol_mass(c)%conc(k,j,i)
                          ENDDO
-                         local_pf(i,j,k) = temp_bin
+                         tend(k,j,i) = temp_bin
                       ENDDO
                    ENDDO
                 ENDDO
              ELSE
-                IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_BC_av
-                IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_DU_av   
-                IF ( TRIM( variable(3:) ) == 'NH' )   to_be_resorted => s_NH_av   
-                IF ( TRIM( variable(3:) ) == 'NO' )   to_be_resorted => s_NO_av   
-                IF ( TRIM( variable(3:) ) == 'OC' )   to_be_resorted => s_OC_av   
-                IF ( TRIM( variable(3:) ) == 'SO4' )  to_be_resorted => s_SO4_av   
-                IF ( TRIM( variable(3:) ) == 'SS' )   to_be_resorted => s_SS_av 
+                tend = 0.0_wp
+             ENDIF
+             IF ( .NOT. mask_surface(mid) )  THEN   
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
-                      DO  k = 1, mask_size_l(mid,3)                   
-                         local_pf(i,j,k) = to_be_resorted(mask_k(mid,k),       &
-                                                    mask_j(mid,j),mask_i(mid,i))
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), &
+                                                 mask_i(mid,i) )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE      
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),&
+                                                                  mask_i(mid,i),&
+                                                                  grid )
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k),&
+                                                      nzt+1 ),&
+                                                 mask_j(mid,j), mask_i(mid,i) )
                       ENDDO
                    ENDDO
                 ENDDO
              ENDIF
+             resorted = .TRUE.
+          ELSE
+             IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_BC_av
+             IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_DU_av   
+             IF ( TRIM( variable(3:) ) == 'NH' )   to_be_resorted => s_NH_av   
+             IF ( TRIM( variable(3:) ) == 'NO' )   to_be_resorted => s_NO_av   
+             IF ( TRIM( variable(3:) ) == 'OC' )   to_be_resorted => s_OC_av   
+             IF ( TRIM( variable(3:) ) == 'SO4' )  to_be_resorted => s_SO4_av   
+             IF ( TRIM( variable(3:) ) == 'SS' )   to_be_resorted => s_SS_av
           ENDIF
        
@@ -11661,4 +10295,38 @@
     END SELECT
     
+    
+    IF ( .NOT. resorted )  THEN
+       IF ( .NOT. mask_surface(mid) )  THEN
+!
+!--       Default masked output    
+          DO  i = 1, mask_size_l(mid,1)
+             DO  j = 1, mask_size_l(mid,2)
+                DO  k = 1, mask_size_l(mid,3)
+                   local_pf(i,j,k) = to_be_resorted( mask_k(mid,k),            &
+                                                     mask_j(mid,j),mask_i(mid,i) )
+                ENDDO
+             ENDDO
+          ENDDO
+       ELSE
+!
+!--       Terrain-following masked output      
+          DO  i = 1, mask_size_l(mid,1)
+             DO  j = 1, mask_size_l(mid,2)
+!
+!--             Get k index of highest horizontal surface
+                topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),     &
+                                                            mask_i(mid,i), grid )
+!
+!--             Save output array
+                DO  k = 1, mask_size_l(mid,3)
+                   local_pf(i,j,k) = to_be_resorted( MIN( topo_top_ind+mask_k(mid,k),&
+                                                          nzt+1 ),             &
+                                                     mask_j(mid,j), mask_i(mid,i) )
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+    
  END SUBROUTINE salsa_data_output_mask
  
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3579)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -122,5 +123,5 @@
 ! - Workaround for sum-up usm arrays in case of restart runs, to avoid program 
 !   crash (MS)
-! 
+!
 ! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
@@ -263,5 +264,5 @@
     USE control_parameters,                                                    &
         ONLY:  air_chemistry, average_count_3d, biometeorology, doav, doav_n,  &
-               land_surface, ocean_mode, rho_surface, urban_surface,           &
+               land_surface, ocean_mode, rho_surface, salsa, urban_surface,    &
                varnamelength
 
@@ -290,5 +291,5 @@
          
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_3d_data_averaging          
+        ONLY:  salsa_3d_data_averaging          
 
     USE surface_mod,                                                           &
@@ -552,7 +553,5 @@
                    CALL radiation_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
-                
-!
-!--             SALSA quantity
+
                 IF ( salsa )  THEN
                    CALL salsa_3d_data_averaging( 'allocate', doav(ii) )
@@ -1202,7 +1201,5 @@
                 CALL radiation_3d_data_averaging( 'sum', doav(ii) )
              ENDIF
-             
-!
-!--          SALSA quantity
+
              IF ( salsa )  THEN
                 CALL salsa_3d_data_averaging( 'sum', doav(ii) )
Index: palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- palm/trunk/SOURCE/swap_timelevel.f90	(revision 3579)
+++ palm/trunk/SOURCE/swap_timelevel.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -162,5 +163,5 @@
         ONLY:  air_chemistry, humidity, land_surface, neutral, ocean_mode,     &
                passive_scalar, simulated_time, timestep_count, urban_surface,  &
-               time_since_reference_point
+               time_since_reference_point, salsa
 
     USE gust_mod,                                                              &
@@ -182,5 +183,5 @@
          
     USE salsa_mod,                                                             &
-        ONLY:  salsa, salsa_swap_timelevel, skip_time_do_salsa             
+        ONLY:  salsa_swap_timelevel, skip_time_do_salsa             
 
     USE turbulence_closure_mod,                                                &
@@ -295,5 +296,5 @@
     IF ( salsa  .AND.  simulated_time >= time_since_reference_point )  THEN
        CALL salsa_swap_timelevel(  MOD( timestep_count, 2 ) )
-    ENDIF  
+    ENDIF     
 
     CALL tcm_swap_timelevel( MOD( timestep_count, 2) )
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3579)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3582)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! Former revisions:
@@ -460,5 +461,6 @@
                ocean_mode, passive_scalar, pt_reference,                       &
                pt_slope_offset, random_heatflux, rans_mode,                    &
-               rans_tke_e, run_coupled, simulated_time, simulated_time_chr,    &
+               rans_tke_e, run_coupled, salsa,                                 &
+               simulated_time, simulated_time_chr,                             &
                skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
                skip_time_do3d, skip_time_domask, skip_time_dopr,               &
@@ -475,5 +477,5 @@
                use_single_reference_value, u_gtrans, v_gtrans,                 &
                virtual_flight, virtual_measurement, wind_turbine,              &
-               ws_scheme_mom, ws_scheme_sca
+               ws_scheme_mom, ws_scheme_sca 
 
     USE cpulog,                                                                &
@@ -542,4 +544,9 @@
               radiation_interaction, radiation_interactions,                   &
               skip_time_do_radiation, time_radiation
+         
+    USE salsa_mod,                                                             &
+        ONLY: aerosol_number, aerosol_mass, nbins, ncc_tot, ngast,             &
+              salsa_boundary_conds, salsa_gas, salsa_gases_from_chem,          &
+              skip_time_do_salsa                      
 
     USE spectra_mod,                                                           &
@@ -574,10 +581,5 @@
         ONLY:  dt_stg_call, dt_stg_adjust, parametrize_inflow_turbulence,      &
                stg_adjust, stg_main, time_stg_adjust, time_stg_call,           &
-               use_syn_turb_gen
-         
-    USE salsa_mod,                                                             &
-        ONLY: aerosol_number, aerosol_mass, nbins, ncc_tot, ngast, salsa,      &
-              salsa_boundary_conds, salsa_gas, salsa_gases_from_chem,          &
-              skip_time_do_salsa           
+               use_syn_turb_gen   
 
     USE user_actions_mod,                                                      &
@@ -916,5 +918,4 @@
              CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
           ENDIF         
-
           CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
 
Index: palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3579)
+++ palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3582)
@@ -20,4 +20,6 @@
 ! Current revisions:
 ! -----------------
+! Move the control parameter "salsa" from salsa_mod to control_parameters 
+! (M. Kurppa)
 ! 
 ! 
@@ -945,5 +947,5 @@
               
        USE salsa_mod,                                                          &
-           ONLY:  salsa, salsa_wrd_local          
+           ONLY:  salsa_wrd_local          
 
        USE surface_mod,                                                        &
@@ -1273,7 +1275,7 @@
        IF ( land_surface )         CALL lsm_wrd_local
        IF ( ocean_mode )           CALL ocean_wrd_local
-       CALL surface_wrd_local
        IF ( radiation )            CALL radiation_wrd_local
        IF ( salsa )                CALL salsa_wrd_local
+       CALL surface_wrd_local
        IF ( urban_surface )        CALL usm_wrd_local
 
