Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4415)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4416)
@@ -21,5 +21,11 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bug fixes and reformatting for the restart data and averaged data output
+! - add missing arrays (averaged data output) in salsa_wrd_local and
+!   salsa_rrd_local
+! - set write_binary_salsa and read_restart_data_salsa to .T. by default
+! - restructure the average arrays for gases and total mass concentrations of
+!   chemical components: set to 4d arrays instead of separate arrays
+! - add allocation checks for averaged data output arrays
 ! 
 ! Former revisions:
@@ -499,8 +505,8 @@
     LOGICAL ::  nesting_offline_salsa   = .TRUE.   !< Apply offline nesting for salsa
     LOGICAL ::  no_insoluble            = .FALSE.  !< Exclude insoluble chemical components
-    LOGICAL ::  read_restart_data_salsa = .FALSE.  !< Read restart data for salsa
+    LOGICAL ::  read_restart_data_salsa = .TRUE.   !< Read restart data for salsa
     LOGICAL ::  salsa_gases_from_chem   = .FALSE.  !< Transfer the gaseous components to SALSA
     LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Include van der Waals and viscous forces in coagulation
-    LOGICAL ::  write_binary_salsa      = .FALSE.  !< read binary for salsa
+    LOGICAL ::  write_binary_salsa      = .TRUE.   !< read binary for salsa
 !
 !-- Process switches: nl* is read from the NAMELIST and is NOT changed.
@@ -827,11 +833,4 @@
 !-- Data output arrays:
 !
-!-- Gases:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_h2so4_av  !< H2SO4
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_hno3_av   !< HNO3
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_nh3_av    !< NH3
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_ocnv_av   !< non-volatile OC
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_ocsv_av   !< semi-volatile OC
-!
 !-- Integrated:
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ldsa_av  !< lung-deposited surface area
@@ -842,17 +841,14 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  pm10_av  !< PM10
 !
-!-- In the particle phase:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_bc_av   !< black carbon
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_du_av   !< dust
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_h2o_av  !< liquid water
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_nh_av   !< ammonia
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_no_av   !< nitrates
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_oc_av   !< org. carbon
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_so4_av  !< sulphates
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_ss_av   !< sea salt
-!
 !-- Bin specific mass and number concentrations:
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  mbins_av  !< bin mas
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  nbins_av  !< bin number
+!
+!-- Gases:
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  salsa_gases_av  !< gases
+!
+!-- In the particle phase:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_h2o_av  !< liquid water
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  s_mass_av  !< mass components
 
 !
@@ -1217,5 +1213,9 @@
        init_gases_type = 1
     ENDIF
-
+!
+!-- If the run is not a restart run, set read_restart_data to .FALSE.
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+       read_restart_data_salsa = .FALSE.
+    ENDIF
 
  END SUBROUTINE salsa_check_parameters
@@ -2907,4 +2907,12 @@
        SELECT CASE ( restart_string(1:length) )
 
+          CASE ( 'aerosol_mass' )
+             DO  ic = 1, ncomponents_mass * nbins_aerosol
+                IF ( k == 1 )  READ ( 13 ) tmp_3d
+                aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             ENDDO
+             found = .TRUE.
+
           CASE ( 'aerosol_number' )
              DO  ib = 1, nbins_aerosol
@@ -2912,22 +2920,109 @@
                 aerosol_number(ib)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
                                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                found = .TRUE.
              ENDDO
-
-          CASE ( 'aerosol_mass' )
-             DO  ic = 1, ncomponents_mass * nbins_aerosol
+             found = .TRUE.
+
+          CASE( 'salsa_gases_av' )
+             IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
+                ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
+             ENDIF
+             DO  ig = 1, ngases_salsa
                 IF ( k == 1 )  READ ( 13 ) tmp_3d
-                aerosol_mass(ic)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
+                salsa_gases_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ig) =                     &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             ENDDO
+             found = .TRUE.
+
+          CASE ( 'ldsa_av' )
+             IF ( .NOT. ALLOCATED( ldsa_av ) )  ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             ldsa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 'mbins_av' )
+             IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
+                ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
+             ENDIF
+             DO  ib = 1, nbins_aerosol
+                IF ( k == 1 )  READ ( 13 ) tmp_3d
+                mbins_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ib) =                           &
                                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
                 found = .TRUE.
              ENDDO
 
-          CASE ( 'salsa_gas' )
-             DO  ig = 1, ngases_salsa
+          CASE ( 'nbins_av' )
+             IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
+                ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
+             ENDIF
+             DO  ib = 1, nbins_aerosol
                 IF ( k == 1 )  READ ( 13 ) tmp_3d
-                salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                    &
+                nbins_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ib) =                           &
                                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
                 found = .TRUE.
              ENDDO
+
+          CASE ( 'ntot_av' )
+             IF ( .NOT. ALLOCATED( ntot_av ) )  ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             ntot_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 'nufp_av' )
+             IF ( .NOT. ALLOCATED( nufp_av ) )  ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             nufp_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 'pm01_av' )
+             IF ( .NOT. ALLOCATED( pm01_av ) )  ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             pm01_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 'pm25_av' )
+             IF ( .NOT. ALLOCATED( pm25_av ) )  ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             pm25_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 'pm10_av' )
+             IF ( .NOT. ALLOCATED( pm10_av ) )  ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             pm10_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                  &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 's_mass_av' )
+             IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
+                ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass) )
+             ENDIF
+             DO  ic = 1, ncomponents_mass
+                IF ( k == 1 )  READ ( 13 ) tmp_3d
+                s_mass_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,ic) =                          &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             ENDDO
+             found = .TRUE.
+
+          CASE ( 's_h2o_av' )
+             IF ( .NOT. ALLOCATED( s_h2o_av ) )  ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             IF ( k == 1 )  READ ( 13 ) tmp_3d
+             s_h2o_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                 &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+             found = .TRUE.
+
+          CASE ( 'salsa_gas' )
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   IF ( k == 1 )  READ ( 13 ) tmp_3d
+                   salsa_gas(ig)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                ENDDO
+                found = .TRUE.
+             ENDIF
 
           CASE DEFAULT
@@ -2960,4 +3055,9 @@
     IF ( write_binary  .AND.  write_binary_salsa )  THEN
 
+       CALL wrd_write_string( 'aerosol_mass' )
+       DO  ic = 1, nbins_aerosol * ncomponents_mass
+          WRITE ( 14 )  aerosol_mass(ic)%conc
+       ENDDO
+
        CALL wrd_write_string( 'aerosol_number' )
        DO  ib = 1, nbins_aerosol
@@ -2965,13 +3065,76 @@
        ENDDO
 
-       CALL wrd_write_string( 'aerosol_mass' )
-       DO  ic = 1, nbins_aerosol * ncomponents_mass
-          WRITE ( 14 )  aerosol_mass(ic)%conc
-       ENDDO
-
-       CALL wrd_write_string( 'salsa_gas' )
-       DO  ig = 1, ngases_salsa
-          WRITE ( 14 )  salsa_gas(ig)%conc
-       ENDDO
+       IF (  .NOT. salsa_gases_from_chem )  THEN
+
+          IF ( ALLOCATED( salsa_gases_av ) )  THEN
+             CALL wrd_write_string( 'salsa_gases_av' )
+             DO  ig = 1, ngases_salsa
+                WRITE ( 14 )  salsa_gases_av(:,:,:,ig)
+             ENDDO
+          ENDIF
+       ENDIF
+
+       IF ( ALLOCATED( ldsa_av ) )  THEN
+       CALL wrd_write_string( 'ldsa_av' )
+       WRITE ( 14 )  ldsa_av
+       ENDIF
+
+       IF ( ALLOCATED( mbins_av ) )  THEN
+          CALL wrd_write_string( 'mbins_av' )
+          DO  ib = 1, nbins_aerosol
+             WRITE ( 14 )  mbins_av(:,:,:,ib)
+          ENDDO
+       ENDIF
+
+       IF ( ALLOCATED( nbins_av ) )  THEN
+          CALL wrd_write_string( 'nbins_av' )
+          DO  ib = 1, nbins_aerosol
+             WRITE ( 14 )  nbins_av(:,:,:,ib)
+          ENDDO
+       ENDIF
+
+       IF ( ALLOCATED( ldsa_av ) )  THEN
+          CALL wrd_write_string( 'ntot_av' )
+          WRITE ( 14 )  ntot_av
+       ENDIF
+
+       IF ( ALLOCATED( nufp_av ) )  THEN
+          CALL wrd_write_string( 'nufp_av' )
+          WRITE ( 14 )  nufp_av
+       ENDIF
+
+       IF ( ALLOCATED( pm01_av ) )  THEN
+          CALL wrd_write_string( 'pm01_av' )
+          WRITE ( 14 )  pm01_av
+       ENDIF
+
+       IF ( ALLOCATED( pm25_av ) )  THEN
+          CALL wrd_write_string( 'pm25_av' )
+          WRITE ( 14 )  pm25_av
+       ENDIF
+
+       IF ( ALLOCATED( pm10_av ) )  THEN
+          CALL wrd_write_string( 'pm10_av' )
+          WRITE ( 14 )  pm10_av
+       ENDIF
+
+       IF ( ALLOCATED( s_mass_av ) )  THEN
+          CALL wrd_write_string( 's_mass_av' )
+          DO  ic = 1, ncomponents_mass
+             WRITE ( 14 )  s_mass_av(:,:,:,ic)
+          ENDDO
+       ENDIF
+
+       IF ( ALLOCATED( s_h2o_av ) )  THEN
+          CALL wrd_write_string( 's_h2o_av' )
+          WRITE ( 14 )  s_h2o_av
+       ENDIF
+
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          CALL wrd_write_string( 'salsa_gas' )
+          DO  ig = 1, ngases_salsa
+             WRITE ( 14 )  salsa_gas(ig)%conc
+          ENDDO
+       ENDIF
 
     ENDIF
@@ -3378,22 +3541,22 @@
              ic = ic+1
           ENDDO
-       ENDIF
-       IF ( prunmode == 1 )  THEN
-          nc_h2o = get_index( prtcl,'H2O' )
-          vc = 8
-          str = ( nc_h2o-1 ) * nbins_aerosol + 1
-          endi = nc_h2o * nbins_aerosol
-          ic = 1
-          DO ss = str, endi
-             aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
-                                             aero_old(ic)%volc(vc) ) * arhoh2o )
-             IF ( k == nzb+1 )  THEN
-                aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
-             ELSEIF ( k == nzt  )  THEN
-                aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
-                aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
-             ENDIF
-             ic = ic+1
-          ENDDO
+          IF ( prunmode == 1 )  THEN
+             nc_h2o = get_index( prtcl,'H2O' )
+             vc = 8
+             str = ( nc_h2o-1 ) * nbins_aerosol + 1
+             endi = nc_h2o * nbins_aerosol
+             ic = 1
+             DO ss = str, endi
+                aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), ( lo_aero(ic)%volc(vc) - &
+                                                aero_old(ic)%volc(vc) ) * arhoh2o )
+                IF ( k == nzb+1 )  THEN
+                   aerosol_mass(ss)%init(k-1) = aerosol_mass(ss)%init(k)
+                ELSEIF ( k == nzt  )  THEN
+                   aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k)
+                   aerosol_mass(ss)%conc(k+1,j,i) = aerosol_mass(ss)%init(k)
+                ENDIF
+                ic = ic+1
+             ENDDO
+          ENDIF
        ENDIF
 !
@@ -9833,4 +9996,10 @@
        ENDIF
 
+    ELSEIF ( var(7:11) == 's_H2O' )  THEN
+       IF ( .NOT. advect_particle_water )  THEN
+          message_string = 'to output s_H2O/s_H2O_av requires that advect_particle_water = .T.'
+          CALL message( 'check_parameters', 'PA0707', 1, 2, 0, 6, 0 )
+       ENDIF
+
     ELSE
        SELECT CASE ( TRIM( var(7:) ) )
@@ -10169,33 +10338,9 @@
           SELECT CASE ( TRIM( variable(7:) ) )
 
-             CASE ( 'g_H2SO4' )
-                IF ( .NOT. ALLOCATED( g_h2so4_av ) )  THEN
-                   ALLOCATE( g_h2so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
+                IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
+                   ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
                 ENDIF
-                g_h2so4_av = 0.0_wp
-
-             CASE ( 'g_HNO3' )
-                IF ( .NOT. ALLOCATED( g_hno3_av ) )  THEN
-                   ALLOCATE( g_hno3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                g_hno3_av = 0.0_wp
-
-             CASE ( 'g_NH3' )
-                IF ( .NOT. ALLOCATED( g_nh3_av ) )  THEN
-                   ALLOCATE( g_nh3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                g_nh3_av = 0.0_wp
-
-             CASE ( 'g_OCNV' )
-                IF ( .NOT. ALLOCATED( g_ocnv_av ) )  THEN
-                   ALLOCATE( g_ocnv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                g_ocnv_av = 0.0_wp
-
-             CASE ( 'g_OCSV' )
-                IF ( .NOT. ALLOCATED( g_ocsv_av ) )  THEN
-                   ALLOCATE( g_ocsv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                g_ocsv_av = 0.0_wp
+                salsa_gases_av = 0.0_wp
 
              CASE ( 'LDSA' )
@@ -10235,15 +10380,9 @@
                 pm10_av = 0.0_wp
 
-             CASE ( 's_BC' )
-                IF ( .NOT. ALLOCATED( s_bc_av ) )  THEN
-                   ALLOCATE( s_bc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+                IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
+                   ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass+1) )
                 ENDIF
-                s_bc_av = 0.0_wp
-
-             CASE ( 's_DU' )
-                IF ( .NOT. ALLOCATED( s_du_av ) )  THEN
-                   ALLOCATE( s_du_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_du_av = 0.0_wp
+                s_mass_av = 0.0_wp
 
              CASE ( 's_H2O' )
@@ -10253,34 +10392,4 @@
                 s_h2o_av = 0.0_wp
 
-             CASE ( 's_NH' )
-                IF ( .NOT. ALLOCATED( s_nh_av ) )  THEN
-                   ALLOCATE( s_nh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_nh_av = 0.0_wp
-
-             CASE ( 's_NO' )
-                IF ( .NOT. ALLOCATED( s_no_av ) )  THEN
-                   ALLOCATE( s_no_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_no_av = 0.0_wp
-
-             CASE ( 's_OC' )
-                IF ( .NOT. ALLOCATED( s_oc_av ) )  THEN
-                   ALLOCATE( s_oc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_oc_av = 0.0_wp
-
-             CASE ( 's_SO4' )
-                IF ( .NOT. ALLOCATED( s_so4_av ) )  THEN
-                   ALLOCATE( s_so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_so4_av = 0.0_wp
-
-             CASE ( 's_SS' )
-                IF ( .NOT. ALLOCATED( s_ss_av ) )  THEN
-                   ALLOCATE( s_ss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                s_ss_av = 0.0_wp
-
              CASE DEFAULT
                 CONTINUE
@@ -10293,190 +10402,120 @@
 
        IF ( variable(7:11) ==  'N_bin' )  THEN
-          READ( variable(12:),* ) char_to_int
-          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
-             ib = char_to_int
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
-          READ( variable(12:),* ) char_to_int
-          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
-             ib = char_to_int
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      temp_bin = 0.0_wp
-                      DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
-                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
-                      ENDDO
-                      mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + temp_bin
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       ELSE
-
-          SELECT CASE ( TRIM( variable(7:) ) )
-
-             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
-
-                vari = TRIM( variable(9:) )  ! remove salsa_g_ from beginning
-
-                SELECT CASE( vari )
-
-                   CASE( 'H2SO4' )
-                      found_index = 1
-                      to_be_resorted => g_h2so4_av
-
-                   CASE( 'HNO3' )
-                      found_index = 2
-                      to_be_resorted => g_hno3_av
-
-                   CASE( 'NH3' )
-                      found_index = 3
-                      to_be_resorted => g_nh3_av
-
-                   CASE( 'OCNV' )
-                      found_index = 4
-                      to_be_resorted => g_ocnv_av
-
-                   CASE( 'OCSV' )
-                      found_index = 5
-                      to_be_resorted => g_ocsv_av
-
-                END SELECT
-
+          IF ( ALLOCATED( nbins_av ) )  THEN
+             READ( variable(12:),* ) char_to_int
+             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
+                ib = char_to_int
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                           &
-                                                 salsa_gas(found_index)%conc(k,j,i)
+                         nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i)
                       ENDDO
                    ENDDO
                 ENDDO
-
-             CASE ( 'LDSA' )
+             ENDIF
+          ENDIF
+
+       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
+          IF ( ALLOCATED( mbins_av ) )  THEN
+             READ( variable(12:),* ) char_to_int
+             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
+                ib = char_to_int
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          temp_bin = 0.0_wp
-                         DO  ib = 1, nbins_aerosol
-   !
-   !--                      Diameter in micrometres
-                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
-   !
-   !--                      Deposition factor: alveolar (use ra_dry)
-                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
-                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
-                                   1.362_wp )**2 ) )
-   !
-   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
-                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
-                                       aerosol_number(ib)%conc(k,j,i)
+                         DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
+                            temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
                          ENDDO
-                         ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin
+                         mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + temp_bin
                       ENDDO
                    ENDDO
                 ENDDO
-
-             CASE ( 'N_UFP' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         temp_bin = 0.0_wp
-                         DO  ib = 1, nbins_aerosol
-                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
-                               temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
-                            ENDIF
-                         ENDDO
-                         nufp_av(k,j,i) = nufp_av(k,j,i) + temp_bin
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'Ntot' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         DO  ib = 1, nbins_aerosol
-                            ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i)
+             ENDIF
+          ENDIF
+       ELSE
+
+          SELECT CASE ( TRIM( variable(7:) ) )
+
+             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
+                IF ( ALLOCATED( salsa_gases_av ) )  THEN
+
+                   vari = TRIM( variable(9:) )  ! remove salsa_g_ from beginning
+
+                   SELECT CASE( vari )
+
+                      CASE( 'H2SO4' )
+                         found_index = 1
+                      CASE( 'HNO3' )
+                         found_index = 2
+                      CASE( 'NH3' )
+                         found_index = 3
+                      CASE( 'OCNV' )
+                         found_index = 4
+                      CASE( 'OCSV' )
+                         found_index = 5
+
+                   END SELECT
+
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            salsa_gases_av(k,j,i,found_index) = salsa_gases_av(k,j,i,found_index)  &
+                                                                + salsa_gas(found_index)%conc(k,j,i)
                          ENDDO
                       ENDDO
                    ENDDO
-                ENDDO
-
-             CASE ( 'PM0.1' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         temp_bin = 0.0_wp
-                         DO  ib = 1, nbins_aerosol
-                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
-                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
-                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
-                               ENDDO
-                            ENDIF
-                         ENDDO
-                         pm01_av(k,j,i) = pm01_av(k,j,i) + temp_bin
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'PM2.5' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         temp_bin = 0.0_wp
-                         DO  ib = 1, nbins_aerosol
-                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
-                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
-                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
-                               ENDDO
-                            ENDIF
-                         ENDDO
-                         pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'PM10' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         temp_bin = 0.0_wp
-                         DO  ib = 1, nbins_aerosol
-                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
-                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
-                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
-                               ENDDO
-                            ENDIF
-                         ENDDO
-                         pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
-                IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
-                   found_index = get_index( prtcl, TRIM( variable(9:) ) )
-                   IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
-                   IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
-                   IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
-                   IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
-                   IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
-                   IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
-                   IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
+                ENDIF
+
+             CASE ( 'LDSA' )
+                IF ( ALLOCATED( ldsa_av ) )  THEN
                    DO  i = nxlg, nxrg
                       DO  j = nysg, nyng
                          DO  k = nzb, nzt+1
-                            DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
-                               to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                     &
-                                                       aerosol_mass(ic)%conc(k,j,i)
+                            temp_bin = 0.0_wp
+                            DO  ib = 1, nbins_aerosol
+!
+!--                            Diameter in micrometres
+                               mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
+!
+!--                            Deposition factor: alveolar (use ra_dry)
+                               df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) + &
+                                      2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) &
+                                      - 1.362_wp )**2 ) )
+!
+!--                            Lung-deposited surface area LDSA (units mum2/cm3)
+                               temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *             &
+                                          aerosol_number(ib)%conc(k,j,i)
+                            ENDDO
+                            ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'N_UFP' )
+                IF ( ALLOCATED( nufp_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            temp_bin = 0.0_wp
+                            DO  ib = 1, nbins_aerosol
+                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
+                                  temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
+                               ENDIF
+                            ENDDO
+                            nufp_av(k,j,i) = nufp_av(k,j,i) + temp_bin
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'Ntot' )
+               IF ( ALLOCATED( ntot_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            DO  ib = 1, nbins_aerosol
+                               ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i)
                             ENDDO
                          ENDDO
@@ -10485,153 +10524,237 @@
                 ENDIF
 
+             CASE ( 'PM0.1' )
+                IF ( ALLOCATED( pm01_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            temp_bin = 0.0_wp
+                            DO  ib = 1, nbins_aerosol
+                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
+                                  DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
+                                     temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
+                                  ENDDO
+                               ENDIF
+                            ENDDO
+                            pm01_av(k,j,i) = pm01_av(k,j,i) + temp_bin
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'PM2.5' )
+                IF ( ALLOCATED( pm25_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            temp_bin = 0.0_wp
+                            DO  ib = 1, nbins_aerosol
+                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
+                                  DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
+                                     temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
+                                  ENDDO
+                               ENDIF
+                            ENDDO
+                            pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'PM10' )
+                IF ( ALLOCATED( pm10_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            temp_bin = 0.0_wp
+                            DO  ib = 1, nbins_aerosol
+                               IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
+                                  DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
+                                     temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
+                                  ENDDO
+                               ENDIF
+                            ENDDO
+                            pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+                IF ( ALLOCATED( s_mass_av ) )  THEN
+                   IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
+                      found_index = get_index( prtcl, TRIM( variable(9:) ) )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nzt+1
+                               DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
+                                  s_mass_av(k,j,i,found_index) = s_mass_av(k,j,i,found_index) +    &
+                                                                 aerosol_mass(ic)%conc(k,j,i)
+                               ENDDO
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                   ENDIF
+                ENDIF
+
              CASE ( 's_H2O' )
-                found_index = get_index( prtcl,'H2O' )
-                to_be_resorted => s_h2o_av
+                IF ( ALLOCATED( s_H2O_av ) )  THEN
+                   found_index = get_index( prtcl,'H2O' )
+                   to_be_resorted => s_h2o_av
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
+                               s_h2o_av(k,j,i) = s_h2o_av(k,j,i) + aerosol_mass(ic)%conc(k,j,i)
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE DEFAULT
+                CONTINUE
+
+          END SELECT
+
+       ENDIF
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       IF ( variable(7:11) ==  'N_bin' )  THEN
+          IF ( ALLOCATED( nbins_av ) )  THEN
+             READ( variable(12:),* ) char_to_int
+             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
+                ib = char_to_int
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
-                            s_h2o_av(k,j,i) = s_h2o_av(k,j,i) + aerosol_mass(ic)%conc(k,j,i)
+                         nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
+          IF ( ALLOCATED( mbins_av ) )  THEN
+             READ( variable(12:),* ) char_to_int
+             IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
+                ib = char_to_int
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+       ELSE
+
+          SELECT CASE ( TRIM( variable(7:) ) )
+
+             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
+                IF ( ALLOCATED( salsa_gases_av ) )  THEN
+                   IF ( TRIM( variable(9:) ) == 'H2SO4' )  THEN  ! 9: remove salsa_g_ from beginning
+                      found_index = 1
+                   ELSEIF ( TRIM( variable(9:) ) == 'HNO3' )  THEN
+                      found_index = 2
+                   ELSEIF ( TRIM( variable(9:) ) == 'NH3' )  THEN
+                      found_index = 3
+                   ELSEIF ( TRIM( variable(9:) ) == 'OCNV' )  THEN
+                      found_index = 4
+                   ELSEIF ( TRIM( variable(9:) ) == 'OCSV' )  THEN
+                      found_index = 5
+                   ENDIF
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            salsa_gases_av(k,j,i,found_index) = salsa_gases_av(k,j,i,found_index)  &
+                                                                / REAL( average_count_3d, KIND=wp )
                          ENDDO
                       ENDDO
                    ENDDO
-                ENDDO
-
-             CASE DEFAULT
-                CONTINUE
-
-          END SELECT
-
-       ENDIF
-
-    ELSEIF ( mode == 'average' )  THEN
-
-       IF ( variable(7:11) ==  'N_bin' )  THEN
-          READ( variable(12:),* ) char_to_int
-          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
-             ib = char_to_int
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
-          READ( variable(12:),* ) char_to_int
-          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
-             ib = char_to_int
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       ELSE
-
-          SELECT CASE ( TRIM( variable(7:) ) )
-
-             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
-                IF ( TRIM( variable(9:) ) == 'H2SO4' )  THEN  ! 9: remove salsa_g_ from beginning
-                   found_index = 1
-                   to_be_resorted => g_h2so4_av
-                ELSEIF ( TRIM( variable(9:) ) == 'HNO3' )  THEN
-                   found_index = 2
-                   to_be_resorted => g_hno3_av
-                ELSEIF ( TRIM( variable(9:) ) == 'NH3' )  THEN
-                   found_index = 3
-                   to_be_resorted => g_nh3_av
-                ELSEIF ( TRIM( variable(9:) ) == 'OCNV' )  THEN
-                   found_index = 4
-                   to_be_resorted => g_ocnv_av
-                ELSEIF ( TRIM( variable(9:) ) == 'OCSV' )  THEN
-                   found_index = 5
-                   to_be_resorted => g_ocsv_av
                 ENDIF
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                           &
-                                                 REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
 
              CASE ( 'LDSA' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'N_UFP' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         nufp_av(k,j,i) = nufp_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'Ntot' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-
-             CASE ( 'PM0.1' )
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         pm01_av(k,j,i) = pm01_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'PM2.5' )
-                DO  i = nxlg, nxrg
-                   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 )
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 'PM10' )
-                DO  i = nxlg, nxrg
-                   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 )
-                      ENDDO
-                   ENDDO
-                ENDDO
-
-             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
-                IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
-                   IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
-                   IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
-                   IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
-                   IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
-                   IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
-                   IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
-                   IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av 
+                IF ( ALLOCATED( ldsa_av ) )  THEN
                    DO  i = nxlg, nxrg
                       DO  j = nysg, nyng
                          DO  k = nzb, nzt+1
-                            to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                        &
-                                                    REAL( average_count_3d, KIND=wp )
+                            ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                          ENDDO
                       ENDDO
                    ENDDO
+                ENDIF
+
+             CASE ( 'N_UFP' )
+                IF ( ALLOCATED( nufp_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            nufp_av(k,j,i) = nufp_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'Ntot' )
+                IF ( ALLOCATED( ntot_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'PM0.1' )
+                IF ( ALLOCATED( pm01_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            pm01_av(k,j,i) = pm01_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'PM2.5' )
+                IF ( ALLOCATED( pm25_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      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 )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 'PM10' )
+                IF ( ALLOCATED( pm10_av ) )  THEN
+                   DO  i = nxlg, nxrg
+                      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 )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+
+             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
+                IF ( ALLOCATED( s_mass_av ) )  THEN
+                   IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
+                      found_index = get_index( prtcl, TRIM( variable(9:) ) )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nzt+1
+                               s_mass_av(k,j,i,found_index) = s_mass_av(k,j,i,found_index) /       &
+                                                              REAL( average_count_3d, KIND=wp )
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                   ENDIF
                 ENDIF
 
@@ -10696,6 +10819,4 @@
 
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< output
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted           !< pointer
 !
 !-- Next statement is to avoid compiler warning about unused variable. May be removed in future.
@@ -10721,4 +10842,8 @@
              ENDDO
           ELSE
+             IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
+                ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
+                nbins_av = REAL( fill_value, KIND = wp )
+             ENDIF
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -10753,4 +10878,8 @@
              ENDDO
           ELSE
+             IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
+                ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
+                mbins_av = REAL( fill_value, KIND = wp )
+             ENDIF
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -10771,10 +10900,10 @@
           CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
              vari = TRIM( variable( 9:LEN( TRIM( variable ) ) - 3 ) )  ! 9: remove salsa_g_
+             IF ( vari == 'H2SO4')  found_index = 1
+             IF ( vari == 'HNO3')   found_index = 2
+             IF ( vari == 'NH3')    found_index = 3
+             IF ( vari == 'OCNV')   found_index = 4
+             IF ( vari == 'OCSV')   found_index = 5
              IF ( av == 0 )  THEN
-                IF ( vari == 'H2SO4')  found_index = 1
-                IF ( vari == 'HNO3')   found_index = 2
-                IF ( vari == 'NH3')    found_index = 3
-                IF ( vari == 'OCNV')   found_index = 4
-                IF ( vari == 'OCSV')   found_index = 5
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -10787,14 +10916,14 @@
                 ENDDO
              ELSE
-                IF ( vari == 'H2SO4' )  to_be_resorted => g_h2so4_av
-                IF ( vari == 'HNO3' )   to_be_resorted => g_hno3_av
-                IF ( vari == 'NH3' )    to_be_resorted => g_nh3_av
-                IF ( vari == 'OCNV' )   to_be_resorted => g_ocnv_av
-                IF ( vari == 'OCSV' )   to_be_resorted => g_ocsv_av
+                IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
+                   ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
+                   salsa_gases_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
                       DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
-                                               KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                         local_pf(i,j,k) = MERGE( salsa_gases_av(k,j,i,found_index),               &
+                                                  REAL( fill_value, KIND = wp ),                   &
+                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
                       ENDDO
                    ENDDO
@@ -10811,14 +10940,14 @@
                          temp_bin = 0.0_wp
                          DO  ib = 1, nbins_aerosol
-   !
-   !--                      Diameter in micrometres
+!
+!--                         Diameter in micrometres
                             mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp 
-   !
-   !--                      Deposition factor: alveolar
+!
+!--                         Deposition factor: alveolar
                             df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
                                    2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
                                    1.362_wp )**2 ) )
-   !
-   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
+!
+!--                         Lung-deposited surface area LDSA (units mum2/cm3)
                             temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
                                        aerosol_number(ib)%conc(k,j,i)
@@ -10831,4 +10960,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( ldsa_av ) )  THEN
+                   ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ldsa_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -10861,4 +10994,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( nufp_av ) )  THEN
+                   ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   nufp_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -10889,4 +11026,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( ntot_av ) )  THEN
+                   ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ntot_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -10920,4 +11061,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( pm01_av ) )  THEN
+                   ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   pm01_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -10951,4 +11096,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( pm25_av ) )  THEN
+                   ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   pm25_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -10982,4 +11131,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( pm10_av ) )  THEN
+                   ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   pm10_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11012,16 +11165,14 @@
                    ENDDO
                 ELSE
-                   IF ( vari == 'BC' )   to_be_resorted => s_bc_av
-                   IF ( vari == 'DU' )   to_be_resorted => s_du_av
-                   IF ( vari == 'NH' )   to_be_resorted => s_nh_av
-                   IF ( vari == 'NO' )   to_be_resorted => s_no_av
-                   IF ( vari == 'OC' )   to_be_resorted => s_oc_av
-                   IF ( vari == 'SO4' )  to_be_resorted => s_so4_av
-                   IF ( vari == 'SS' )   to_be_resorted => s_ss_av
+                   IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
+                      ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass) )
+                      s_mass_av = REAL( fill_value, KIND = wp )
+                   ENDIF
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
                          DO  k = nzb_do, nzt_do
-                            local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,      &
-                                               KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                            local_pf(i,j,k) = MERGE( s_mass_av(k,j,i,found_index),                 &
+                                                     REAL( fill_value, KIND = wp ),                &
+                                                     BTEST( wall_flags_total_0(k,j,i), 0 ) )
                          ENDDO
                       ENDDO
@@ -11050,10 +11201,14 @@
                 ENDDO
              ELSE
-                to_be_resorted => s_h2o_av
+     !           to_be_resorted => s_h2o_av
+                IF ( .NOT. ALLOCATED( s_h2o_av ) )  THEN
+                   ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   s_h2o_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
                       DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
-                                              KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                         local_pf(i,j,k) = MERGE( s_h2o_av(k,j,i), REAL( fill_value, KIND = wp ),  &
+                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
                       ENDDO
                    ENDDO
@@ -11110,6 +11265,4 @@
 
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< local
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< pointer
 
     found     = .TRUE.
@@ -11131,4 +11284,8 @@
              ENDDO
           ELSE
+             IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
+                ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
+                nbins_av = REAL( fill_value, KIND = wp )
+             ENDIF
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -11161,4 +11318,8 @@
              ENDDO
           ELSE
+             IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
+                ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
+                mbins_av = REAL( fill_value, KIND = wp )
+             ENDIF
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -11176,11 +11337,11 @@
 
           CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
+             IF ( TRIM( variable(7:) ) == 'g_H2SO4')  found_index = 1
+             IF ( TRIM( variable(7:) ) == 'g_HNO3')   found_index = 2
+             IF ( TRIM( variable(7:) ) == 'g_NH3')    found_index = 3
+             IF ( TRIM( variable(7:) ) == 'g_OCNV')   found_index = 4
+             IF ( TRIM( variable(7:) ) == 'g_OCSV')   found_index = 5
+
              IF ( av == 0 )  THEN
-                IF ( TRIM( variable(7:) ) == 'g_H2SO4')  found_index = 1
-                IF ( TRIM( variable(7:) ) == 'g_HNO3')   found_index = 2
-                IF ( TRIM( variable(7:) ) == 'g_NH3')    found_index = 3
-                IF ( TRIM( variable(7:) ) == 'g_OCNV')   found_index = 4
-                IF ( TRIM( variable(7:) ) == 'g_OCSV')   found_index = 5
-
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11193,16 +11354,16 @@
                 ENDDO
              ELSE
-!
-!--             9: remove salsa_g_ from the beginning
-                IF ( TRIM( variable(9:) ) == 'H2SO4' ) to_be_resorted => g_h2so4_av
-                IF ( TRIM( variable(9:) ) == 'HNO3' )  to_be_resorted => g_hno3_av
-                IF ( TRIM( variable(9:) ) == 'NH3' )   to_be_resorted => g_nh3_av
-                IF ( TRIM( variable(9:) ) == 'OCNV' )  to_be_resorted => g_ocnv_av
-                IF ( TRIM( variable(9:) ) == 'OCSV' )  to_be_resorted => g_ocsv_av
+                IF ( .NOT. ALLOCATED( salsa_gases_av ) )  THEN
+                   ALLOCATE( salsa_gases_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
+                   salsa_gases_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
                       DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
-                                               KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                          local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
+!                                                KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                         local_pf(i,j,k) = MERGE( salsa_gases_av(k,j,i,found_index),               &
+                                                  REAL( fill_value, KIND = wp ),                   &
+                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
                       ENDDO
                    ENDDO
@@ -11217,14 +11378,14 @@
                          temp_bin = 0.0_wp
                          DO  ib = 1, nbins_aerosol
-   !
-   !--                      Diameter in micrometres
+!
+!--                         Diameter in micrometres
                             mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
-   !
-   !--                      Deposition factor: alveolar
+!
+!--                         Deposition factor: alveolar
                             df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
                                    2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
                                    1.362_wp )**2 ) )
-   !
-   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
+!
+!--                         Lung-deposited surface area LDSA (units mum2/cm3)
                             temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
                                        aerosol_number(ib)%conc(k,j,i)
@@ -11236,4 +11397,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( ldsa_av ) )  THEN
+                   ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ldsa_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11263,4 +11428,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( nufp_av ) )  THEN
+                   ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   nufp_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11288,4 +11457,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( ntot_av ) )  THEN
+                   ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ntot_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11317,4 +11490,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( pm01_av ) )  THEN
+                   ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   pm01_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11346,4 +11523,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( pm25_av ) )  THEN
+                   ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   pm25_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11375,4 +11556,8 @@
                 ENDDO
              ELSE
+                IF ( .NOT. ALLOCATED( pm10_av ) )  THEN
+                   ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   pm10_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -11402,18 +11587,14 @@
                    ENDDO
                 ELSE
-!
-!--                9: remove salsa_s_ from the beginning
-                   IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
-                   IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
-                   IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
-                   IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
-                   IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
-                   IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
-                   IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
+                   IF ( .NOT. ALLOCATED( s_mass_av ) )  THEN
+                      ALLOCATE( s_mass_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass) )
+                      s_mass_av = REAL( fill_value, KIND = wp )
+                   ENDIF
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
                          DO  k = nzb_do, nzt_do
-                            local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,      &
-                                                     KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                            local_pf(i,j,k) = MERGE( s_mass_av(k,j,i,found_index),                 &
+                                                     REAL( fill_value, KIND = wp ),                &
+                                                     BTEST( wall_flags_total_0(k,j,i), 0 ) )
                          ENDDO
                       ENDDO
@@ -11438,10 +11619,13 @@
                 ENDDO
              ELSE
-                to_be_resorted => s_h2o_av
+                IF ( .NOT. ALLOCATED( s_h2o_av ) )  THEN
+                   ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   s_h2o_av = REAL( fill_value, KIND = wp )
+                ENDIF
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
                       DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
-                                               KIND = wp ), BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                         local_pf(i,j,k) = MERGE( s_h2o_av(k,j,i), REAL( fill_value, KIND = wp ),  &
+                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
                       ENDDO
                    ENDDO
@@ -11622,9 +11806,10 @@
                 IF ( vari == 'g_OCSV')   to_be_resorted => salsa_gas(5)%conc
              ELSE
-                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
+                IF ( vari == 'g_H2SO4') temp_array = salsa_gases_av(:,:,:,1)
+                IF ( vari == 'g_HNO3')  temp_array = salsa_gases_av(:,:,:,2)
+                IF ( vari == 'g_NH3')   temp_array = salsa_gases_av(:,:,:,3)
+                IF ( vari == 'g_OCNV')  temp_array = salsa_gases_av(:,:,:,4)
+                IF ( vari == 'g_OCSV')  temp_array = salsa_gases_av(:,:,:,5)
+                to_be_resorted => temp_array
              ENDIF
 
@@ -11636,14 +11821,14 @@
                          temp_bin = 0.0_wp
                          DO  ib = 1, nbins_aerosol
-   !
-   !--                      Diameter in micrometres
+!
+!--                         Diameter in micrometres
                             mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
-   !
-   !--                      Deposition factor: alveolar
+!
+!--                         Deposition factor: alveolar
                             df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
                                    2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
                                    1.362_wp )**2 ) )
-   !
-   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
+!
+!--                         Lung-deposited surface area LDSA (units mum2/cm3)
                             temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
                                        aerosol_number(ib)%conc(k,j,i)
@@ -11854,5 +12039,5 @@
                       ENDDO
                    ENDDO
-                ENDDO 
+                ENDDO
                 IF ( .NOT. mask_surface(mid) )  THEN
                    DO  i = 1, mask_size_l(mid,1)
@@ -11943,7 +12128,7 @@
 
           CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
-             IF ( av == 0 )  THEN
-                IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN
-                   found_index = get_index( prtcl, TRIM( variable(9:) ) )
+             IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN
+                found_index = get_index( prtcl, TRIM( variable(9:) ) )
+                IF ( av == 0 )  THEN
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
@@ -11957,48 +12142,41 @@
                       ENDDO
                    ENDDO
-                ELSE
-                   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) = tend( 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)
-!
-!--                      Get k index of the highest terraing surface
-                         im = mask_i(mid,i)
-                         jm = mask_j(mid,j)
-                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
-                                                       DIM = 1 ) - 1
-                         DO  k = 1, mask_size_l(mid,3)
-                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
-!
-!--                         Set value if not in building
-                            IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
-                               local_pf(i,j,k) = fill_value
-                            ELSE
-                               local_pf(i,j,k) = tend(kk,jm,im)
-                            ENDIF
+                   ELSE
+                      DO  i = 1, mask_size_l(mid,1)
+                         DO  j = 1, mask_size_l(mid,2)
+   !
+   !--                      Get k index of the highest terraing surface
+                            im = mask_i(mid,i)
+                            jm = mask_j(mid,j)
+                            ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
+                                                          DIM = 1 ) - 1
+                            DO  k = 1, mask_size_l(mid,3)
+                               kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+   !
+   !--                         Set value if not in building
+                               IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
+                                  local_pf(i,j,k) = fill_value
+                               ELSE
+                                  local_pf(i,j,k) = tend(kk,jm,im)
+                               ENDIF
+                            ENDDO
                          ENDDO
                       ENDDO
-                   ENDDO
+                   ENDIF
+                   resorted = .TRUE.
+                ELSE
+                   temp_array = s_mass_av(:,:,:,found_index)
+                   to_be_resorted => temp_array
                 ENDIF
-                resorted = .TRUE.
              ELSE
-!
-!--             9: remove salsa_s_ from the beginning
-                IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
-                IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
-                IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
-                IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
-                IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
-                IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
-                IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
+                local_pf = fill_value
              ENDIF
 
@@ -12074,4 +12252,5 @@
           DO  i = 1, mask_size_l(mid,1)
              DO  j = 1, mask_size_l(mid,2)
+!
 !--             Get k index of the highest terraing surface
                 im = mask_i(mid,i)
@@ -12081,4 +12260,5 @@
                 DO  k = 1, mask_size_l(mid,3)
                    kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
 !--                Set value if not in building
                    IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
