Index: palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3885)
+++ palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3968)
@@ -27,4 +27,10 @@
 ! -----------------
 ! $Id$
+! - in subroutine chem_emissions_match replace all decision structures relating to
+!   mode_emis to emiss_lod
+! - in subroutine chem_check_parameters replace emt%nspec with emt%n_emiss_species
+! - spring cleaning (E.C. Chan)
+! 
+! 3885 2019-04-11 11:29:34Z kanani
 ! Changes related to global restructuring of location messages and introduction 
 ! of additional debug messages
@@ -36,10 +42,10 @@
 ! Removed unused variables from chem_emissions_mod
 ! 
-!3772 2019-02-28 15:51:57Z suehring 
+! 3772 2019-02-28 15:51:57Z suehring 
 ! - In case of parametrized emissions, assure that emissions are only on natural
 !   surfaces (i.e. streets) and not on urban surfaces. 
 ! - some unnecessary if clauses removed
 !
-!3685 2019 -01-21 01:02:11Z knoop
+! 3685 2019 -01-21 01:02:11Z knoop
 ! Some interface calls moved to module_interface + cleanup
 ! 
@@ -172,7 +178,5 @@
     REAL(wp), PARAMETER ::  r_cp    = 0.286_wp               !< R / cp (exponent for potential temperature)
 
-
     SAVE
-
 
 !    
@@ -196,6 +200,4 @@
        MODULE PROCEDURE chem_emissions_setup
     END INTERFACE chem_emissions_setup
-
-
     
     PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
@@ -213,12 +215,14 @@
  SUBROUTINE chem_emissions_check_parameters
 
-
     IMPLICIT NONE
 
     TYPE(chem_emis_att_type) ::  emt
 
-    !
-    !-- Check Emission Species Number equal to number of passed names for the chemistry species: 
-    IF ( SIZE(emt%species_name) /= emt%nspec  )  THEN
+!
+!-- Check if species count matches the number of names
+!-- passed for the chemiscal species
+
+    IF  ( SIZE(emt%species_name) /= emt%n_emiss_species  )  THEN
+!    IF  ( SIZE(emt%species_name) /= emt%nspec  )  THEN
 
        message_string = 'Numbers of input emission species names and number of species'     //      &
@@ -227,398 +231,64 @@
             
     ENDIF
-    
 
  END SUBROUTINE chem_emissions_check_parameters
+
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Matching the chemical species indices. The routine checks what are the indices of the emission input species 
-!> and the corresponding ones of the model species. The routine gives as output a vector containing the number 
-!> of common species: it is important to note that while the model species are distinct, their values could be 
-!> given to a single species in input: for example, in the case of NO2 and NO, values may be passed in input as
-!> NOx values.
+!> Matching the chemical species indices. The routine checks what are the
+!> indices of the emission input species and the corresponding ones of the
+!> model species. The routine gives as output a vector containing the number 
+!> of common species: it is important to note that while the model species
+!> are distinct, their values could be given to a single species in input.
+!> For example, in the case of NO2 and NO, values may be passed in input as
+!> NOX values.
 !------------------------------------------------------------------------------!
+
 SUBROUTINE chem_emissions_match( emt_att,len_index )   
 
 
-    INTEGER(iwp), INTENT(INOUT)             ::  len_index   !< Variable where to store the number of common species between the input dataset and the model species   
-
-    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att     !< Chemistry Emission Array containing information for all the input chemical emission species
-   
-    INTEGER(iwp) ::  ind_mod, ind_inp                       !< Parameters for cycling through chemical model and input species 
-    INTEGER(iwp) ::  nspec_emis_inp                         !< Variable where to store the number of the emission species in input
-    INTEGER(iwp) ::  ind_voc                                !< Indices to check whether a split for voc should be done
-    INTEGER(iwp) ::  ispec                                  !< index for cycle over effective number of emission species
-
-
-    IF ( debug_output )  CALL debug_message( 'chem_emissions_match', 'start' )
-
-    !
-    !-- Number of input emission species.
-    nspec_emis_inp=emt_att%nspec 
-
-    !
-    !-- Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED 
-    SELECT CASE( TRIM( mode_emis ) ) 
-
-       !
-       !-- PRE-PROCESSED mode
-       CASE ( "PRE-PROCESSED" )
+
+    INTEGER(iwp)  ::  ind_inp                    !< Parameters for cycling through chemical input species 
+    INTEGER(iwp)  ::  ind_mod                    !< Parameters for cycling through chemical model species 
+    INTEGER(iwp)  ::  ind_voc                    !< Indices to check whether a split for voc should be done
+    INTEGER(iwp)  ::  ispec                      !< index for cycle over effective number of emission species
+    INTEGER(iwp)  ::  nspec_emis_inp             !< Variable where to store # of emission species in input
+
+    INTEGER(iwp), INTENT(INOUT)  ::  len_index   !< number of common species between input dataset & model species
+
+    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att     !< Chemistry Emission Array (decl. netcdf_data_input.f90)
+
+
+    IF  ( debug_output )  CALL debug_message( 'chem_emissions_match', 'start' )
+
+!
+!-- Number of input emission species
+
+    nspec_emis_inp = emt_att%n_emiss_species
+!    nspec_emis_inp=emt_att%nspec 
+
+!
+!-- Check the emission LOD: 0 (PARAMETERIZED), 1 (DEFAULT), 2 (PRE-PROCESSED)
+!
+    SELECT CASE (emiss_lod)
+
+!
+!-- LOD 0 (PARAMETERIZED mode)
+
+       CASE (0)
 
           len_index = 0
-          len_index_voc = 0
-
-          IF ( nvar > 0 .AND. (nspec_emis_inp > 0) )  THEN
-             !
-             !-- Cycle over model species
-             DO ind_mod = 1,  nvar 
-                !
-                !-- Cycle over input species  
-                DO ind_inp = 1, nspec_emis_inp
-
-                   !
-                   !-- Check for VOC Species  
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN        
-                      DO ind_voc = 1, emt_att%nvoc
-             
-                         IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                            len_index = len_index + 1
-                            len_index_voc = len_index_voc + 1
-                         ENDIF
-                      END DO
-                   ENDIF
-                   !
-                   !-- Other Species
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                      len_index = len_index + 1
-                   ENDIF
-                ENDDO
-             ENDDO
-
-             !
-             !-- Allocate array for storing the indices of the matched species
-             IF ( len_index > 0 )  THEN
- 
-                ALLOCATE( match_spec_input(len_index) ) 
- 
-                ALLOCATE( match_spec_model(len_index) )
-
-                IF ( len_index_voc > 0 )  THEN
-                   !
-                   !-- contains indices of the VOC model species
-                   ALLOCATE( match_spec_voc_model(len_index_voc) )
-                   !
-                   !-- contains the indices of different values of VOC composition of input variable VOC_composition
-                   ALLOCATE( match_spec_voc_input(len_index_voc) )
-
-                ENDIF
-
-                !
-                !-- pass the species indices to declared arrays
-                len_index = 0
-
-                !
-                !-- Cycle over model species
-                DO ind_mod = 1, nvar 
-                   !
-                   !-- Cycle over Input species  
-                   DO ind_inp = 1, nspec_emis_inp
-                      !
-                      !-- VOCs
-                      IF ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND.    &
-                           ALLOCATED(match_spec_voc_input) )  THEN
-                         
-                         DO ind_voc= 1, emt_att%nvoc
-                            IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                               len_index = len_index + 1
-                               len_index_voc = len_index_voc + 1
-                        
-                               match_spec_input(len_index) = ind_inp
-                               match_spec_model(len_index) = ind_mod
-
-                               match_spec_voc_input(len_index_voc) = ind_voc
-                               match_spec_voc_model(len_index_voc) = ind_mod                         
-                            ENDIF
-                         END DO
-                      ENDIF
-
-                      !
-                      !-- Other Species
-                      IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                         len_index = len_index + 1
-                         match_spec_input(len_index) = ind_inp
-                         match_spec_model(len_index) = ind_mod
-                      ENDIF
-                   END DO
-                END DO
-
-             ELSE
-                !
-                !-- in case there are no species matching
-                message_string = 'Non of given emission species'            //         &
-                                 ' matches'                                //          &
-                                 ' model chemical species:'                //          &
-                                 ' Emission routine is not called'  
-                CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
-             ENDIF
- 
-          ELSE
-
-             !
-             !-- either spc_names is zero or nspec_emis_inp is not allocated
-             message_string = 'Array of Emission species not allocated:'             //          &
-                              ' Either no emission species are provided as input or'  //         &
-                              ' no chemical species are used by PALM:'                //         &
-                              ' Emission routine is not called'                  
-             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) 
-
-          ENDIF
-
-       !
-       !-- DEFAULT mode
-       CASE ("DEFAULT")
-
-          len_index = 0          !<  index for TOTAL number of species   
-          len_index_voc = 0      !<  index for TOTAL number of VOCs
-          len_index_pm = 3       !<  index for TOTAL number of PMs: PM1, PM2.5, PM10.
- 
-          IF ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
-
-             !
-             !-- Cycle over model species
-             DO ind_mod = 1, nvar
-                !
-                !-- Cycle over input species
-                DO ind_inp = 1, nspec_emis_inp
-
-                   !
-                   !-- Check for VOC Species  
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
-                      DO ind_voc= 1, emt_att%nvoc
-                        
-                         IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                            len_index = len_index + 1
-                            len_index_voc = len_index_voc + 1
-                         ENDIF
-                         
-                      END DO
-                   ENDIF
-
-                   !
-                   !-- PMs: There is one input species name for all PM 
-                   !-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
-                      !
-                      !-- PM1
-                      IF ( TRIM( spc_names(ind_mod) ) == "PM1" )  THEN
-                         len_index = len_index + 1
-                      !
-                      !-- PM2.5
-                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
-                         len_index = len_index + 1
-                      !
-                      !-- PM10
-                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
-                         len_index = len_index + 1
-                      ENDIF
-                   ENDIF
-
-                   !
-                   !-- NOx: NO2 and NO    
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" )  THEN 
-                      !
-                      !-- NO
-                      IF ( TRIM( spc_names(ind_mod) ) == "NO" )  THEN
-                         len_index = len_index + 1
-                      !
-                      !-- NO2
-                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
-                         len_index = len_index + 1
-                      ENDIF
-                   ENDIF
-
-                   !
-                   !-- SOx: SO2 and SO4
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" )  THEN
-                      !
-                      !-- SO2
-                      IF ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
-                         len_index = len_index + 1
-                      !
-                      !-- SO4
-                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
-                         len_index = len_index + 1
-                      ENDIF
-                   ENDIF
-
-                   !
-                   !-- Other Species
-                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                      len_index = len_index + 1
-                   ENDIF
-                END DO
-             END DO
-
-
-             !
-             !-- Allocate arrays
-             IF ( len_index > 0 )  THEN
-
-                ALLOCATE( match_spec_input(len_index) )
-                ALLOCATE( match_spec_model(len_index) )
-
-                IF ( len_index_voc > 0 )  THEN
-                   !
-                   !-- Contains indices of the VOC model species
-                   ALLOCATE( match_spec_voc_model(len_index_voc) )
-                   !
-                   !-- Contains the indices of different values of VOC composition
-                   !-- of input variable VOC_composition
-                   ALLOCATE( match_spec_voc_input(len_index_voc) )                                                  
-                ENDIF
-
-                !
-                !-- Pass the species indices to declared arrays
-                len_index = 0
-                len_index_voc = 0
-                
-                DO ind_mod = 1, nvar
-                   DO ind_inp = 1, nspec_emis_inp 
-
-                      !
-                      !-- VOCs
-                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND.   &
-                         ALLOCATED(match_spec_voc_input) )  THEN      
-                         DO ind_voc= 1, emt_att%nvoc
-                            IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                               len_index = len_index + 1
-                               len_index_voc = len_index_voc + 1
-                        
-                               match_spec_input(len_index) = ind_inp
-                               match_spec_model(len_index) = ind_mod
-
-                               match_spec_voc_input(len_index_voc) = ind_voc
-                               match_spec_voc_model(len_index_voc) = ind_mod                         
-                            ENDIF
-                         END DO
-                      ENDIF
-
-                      !
-                      !-- PMs 
-                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
-                         !
-                         !-- PM1
-                         IF ( TRIM( spc_names(ind_mod) ) == "PM1" )  THEN
-                            len_index = len_index + 1
-
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-                         !
-                         !-- PM2.5
-                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
-                            len_index = len_index + 1
-
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-                         !
-                         !-- PM10
-                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
-                            len_index = len_index + 1
-    
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-  
-                         ENDIF
-                      ENDIF
-
-                      !
-                      !-- NOx 
-                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" )  THEN
-                         !
-                         !-- NO
-                         IF ( TRIM( spc_names(ind_mod) ) == "NO" )  THEN
-                            len_index = len_index + 1
-
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-                         !
-                         !-- NO2
-                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
-                            len_index = len_index + 1
-
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-  
-                         ENDIF
-                      ENDIF
-
-                      !
-                      !-- SOx 
-                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" ) THEN
-                         !
-                         !-- SO2
-                         IF ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
-                            len_index = len_index + 1
-
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-  
-                         !
-                         !-- SO4
-                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
-                            len_index = len_index + 1
-
-                            match_spec_input(len_index) = ind_inp
-                            match_spec_model(len_index) = ind_mod
-  
-                         ENDIF
-                      ENDIF
-
-                      !
-                      !-- Other Species
-                      IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
-                         len_index = len_index + 1
-                            
-                         match_spec_input(len_index) = ind_inp
-                         match_spec_model(len_index) = ind_mod
-                      ENDIF
-                   END DO
-                END DO
-
-             ELSE
-
-                message_string = 'Non of given Emission Species'            //         &
-                                 ' matches'                                //          &
-                                 ' model chemical species'                 //          &
-                                 ' Emission routine is not called'         
-                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 
-
-             ENDIF
-
-          ELSE
-
-             message_string = 'Array of Emission species not allocated: '            //          &
-                              ' Either no emission species are provided as input or'  //         &
-                              ' no chemical species are used by PALM:'                //         &
-                              ' Emission routine is not called'                                   
-             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) 
- 
-          ENDIF
- 
-       !
-       !-- PARAMETERIZED mode
-       CASE ("PARAMETERIZED")
-
-          len_index = 0
-
-          IF ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
-
-             !
-             !-- Cycle over model species
+
+          IF  ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
+
+!
+!-- Cycle over model species
+
              DO ind_mod = 1, nvar
                 ind_inp = 1
                 DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )    !< 'novalue' is the default  
-                   IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                   IF  ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                       len_index = len_index + 1
                    ENDIF
@@ -627,20 +297,22 @@
              ENDDO
 
-             IF ( len_index > 0 )  THEN
-
-                !
-                !-- Allocation of Arrays of the matched species
-                ALLOCATE( match_spec_input(len_index) ) 
- 
-                ALLOCATE( match_spec_model(len_index) )
-
-                !
-                !-- Pass the species indices to declared arrays    
+             IF  ( len_index > 0 )  THEN
+
+!
+!-- Allocation of Arrays of the matched species
+
+                ALLOCATE ( match_spec_input(len_index) ) 
+                ALLOCATE ( match_spec_model(len_index) )
+
+!
+!-- Pass species indices to declared arrays
+
                 len_index = 0
 
                 DO ind_mod = 1, nvar  
                    ind_inp = 1
-                   DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )     
-                      IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                   DO WHILE  ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
+                      IF  ( TRIM(surface_csflux_name(ind_inp)) ==          &
+                            TRIM(spc_names(ind_mod))            )  THEN
                          len_index = len_index + 1
                          match_spec_input(len_index) = ind_inp
@@ -651,12 +323,13 @@
                 END DO
 
-                !
-                !-- Check
+!
+!-- Check
+
                 DO ispec = 1, len_index
 
-                   IF ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND.    &
+                   IF  ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND.    &
                         emiss_factor_side(match_spec_input(ispec) ) < 0 )  THEN
 
-                      message_string = 'PARAMETERIZED emissions mode selected:'            //           &
+                      message_string = 'PARAMETERIZED emissions mode selected:'             //          &
                                        ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
                                        ' but values of scaling factors for street types'    //          &
@@ -665,6 +338,8 @@
                                        ' or not provided at all: PLEASE set a finite value' //          &
                                        ' for these parameters in the chemistry namelist'          
-                      CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 ) 
+                      CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
+ 
                    ENDIF
+
                 END DO
 
@@ -672,14 +347,15 @@
              ELSE
                 
-                message_string = 'Non of given Emission Species'            //          &
+                message_string = 'Non of given Emission Species'           //           &
                                  ' matches'                                //           &
                                  ' model chemical species'                 //           &
                                  ' Emission routine is not called'         
-                CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 ) 
+                CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
+ 
              ENDIF
 
           ELSE
      
-             message_string = 'Array of Emission species not allocated: '            //           &
+             message_string = 'Array of Emission species not allocated: '             //          &
                               ' Either no emission species are provided as input or'  //          &
                               ' no chemical species are used by PALM.'                //          &
@@ -687,13 +363,391 @@
              CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) 
  
-          ENDIF             
-
-
-       !
-       !-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
+          ENDIF
+
+!
+!-- LOD 1 (DEFAULT mode)
+
+       CASE (1)
+
+          len_index = 0          ! TOTAL number of species (to be accumulated)  
+          len_index_voc = 0      ! TOTAL number of VOCs (to be accumulated)
+          len_index_pm = 3       ! TOTAL number of PMs: PM1, PM2.5, PM10.
+ 
+          IF  ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
+
+!
+!-- Cycle over model species
+             DO ind_mod = 1, nvar
+
+!
+!-- Cycle over input species
+
+                DO ind_inp = 1, nspec_emis_inp
+
+!
+!-- Check for VOC Species
+
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
+                      DO ind_voc= 1, emt_att%nvoc
+                        
+                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                            len_index = len_index + 1
+                            len_index_voc = len_index_voc + 1
+                         ENDIF
+                         
+                      END DO
+                   ENDIF
+
+!
+!-- PMs: There is one input species name for all PM 
+!-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
+
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
+
+                      IF      ( TRIM( spc_names(ind_mod) ) == "PM1" )   THEN
+                         len_index = len_index + 1
+                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
+                         len_index = len_index + 1
+                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
+                         len_index = len_index + 1
+                      ENDIF
+
+                   ENDIF
+
+!
+!-- NOX: NO2 and NO
+
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN 
+
+                      IF     ( TRIM( spc_names(ind_mod) ) == "NO"  )  THEN
+                         len_index = len_index + 1
+                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
+                         len_index = len_index + 1
+                      ENDIF
+
+                   ENDIF
+
+!
+!-- SOX: SO2 and SO4
+
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" )  THEN
+
+                      IF      ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
+                         len_index = len_index + 1
+                      ELSEIF  ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
+                         len_index = len_index + 1
+                      ENDIF
+
+                   ENDIF
+
+!
+!-- Other Species
+
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) ==             &
+                        TRIM( spc_names(ind_mod) ) )  THEN
+                      len_index = len_index + 1
+                   ENDIF
+
+                END DO  ! ind_inp ...
+
+             END DO  ! ind_mod ...
+
+
+!
+!-- Allocate arrays
+
+             IF  ( len_index > 0 )  THEN
+
+                ALLOCATE ( match_spec_input(len_index) )
+                ALLOCATE ( match_spec_model(len_index) )
+
+                IF  ( len_index_voc > 0 )  THEN
+
+!
+!-- Contains indices of the VOC model species
+
+                   ALLOCATE( match_spec_voc_model(len_index_voc) )
+
+!
+!-- Contains the indices of different values of VOC composition
+!-- of input variable VOC_composition
+
+                   ALLOCATE( match_spec_voc_input(len_index_voc) )
+
+                ENDIF
+
+!
+!-- Pass the species indices to declared arrays
+
+                len_index = 0
+                len_index_voc = 0
+                
+                DO ind_mod = 1, nvar
+                   DO ind_inp = 1, nspec_emis_inp 
+
+!
+!-- VOCs
+
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND.       &
+                           ALLOCATED (match_spec_voc_input) )  THEN
+
+                         DO ind_voc = 1, emt_att%nvoc
+
+                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) ==                  &
+                                 TRIM( spc_names(ind_mod) ) )  THEN
+
+                               len_index     = len_index + 1
+                               len_index_voc = len_index_voc + 1
+                        
+                               match_spec_input(len_index) = ind_inp
+                               match_spec_model(len_index) = ind_mod
+
+                               match_spec_voc_input(len_index_voc) = ind_voc
+                               match_spec_voc_model(len_index_voc) = ind_mod
+
+                            ENDIF
+
+                         END DO
+
+                      ENDIF
+
+!
+!-- PMs
+
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
+
+                         IF      ( TRIM( spc_names(ind_mod) ) == "PM1"  )  THEN
+                            len_index = len_index + 1
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
+                            len_index = len_index + 1
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
+                            len_index = len_index + 1
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+                         ENDIF
+
+                      ENDIF
+
+!
+!-- NOX 
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN
+
+                         IF      ( TRIM( spc_names(ind_mod) ) == "NO"  )  THEN
+                            len_index = len_index + 1
+
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+
+                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
+                            len_index = len_index + 1
+
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+  
+                         ENDIF
+
+                      ENDIF
+
+
+!
+!-- SOX 
+
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
+
+                         IF  ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
+                            len_index = len_index + 1
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+                         ELSEIF  ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
+                            len_index = len_index + 1
+                            match_spec_input(len_index) = ind_inp
+                            match_spec_model(len_index) = ind_mod
+                         ENDIF
+
+                      ENDIF
+
+!
+!-- Other Species
+
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                         len_index = len_index + 1
+                         match_spec_input(len_index) = ind_inp
+                         match_spec_model(len_index) = ind_mod
+                      ENDIF
+
+                   END DO  ! inp_ind
+
+                END DO     ! inp_mod
+
+!
+!-- Error reporting (no matching)
+
+             ELSE
+
+                message_string = 'None of given Emission Species matches'           //   &
+                                 ' model chemical species'                          //   &
+                                 ' Emission routine is not called'         
+                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 
+
+             ENDIF
+
+!
+!-- Error reporting (no species)
+
+          ELSE
+
+             message_string = 'Array of Emission species not allocated: '             // &
+                              ' Either no emission species are provided as input or'  // &
+                              ' no chemical species are used by PALM:'                // &
+                              ' Emission routine is not called'                                   
+             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) 
+ 
+          ENDIF
+
+!
+!-- LOD 2 (PRE-PROCESSED mode)
+
+       CASE (2)
+
+          len_index = 0
+          len_index_voc = 0
+
+          IF  ( nvar > 0 .AND. (nspec_emis_inp > 0) )  THEN
+!
+!-- Cycle over model species
+             DO ind_mod = 1, nvar 
+
+!
+!-- Cycle over input species  
+                DO ind_inp = 1, nspec_emis_inp
+
+!
+!-- Check for VOC Species 
+
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN        
+                      DO ind_voc = 1, emt_att%nvoc
+                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                            len_index     = len_index + 1
+                            len_index_voc = len_index_voc + 1
+                         ENDIF
+                      END DO
+                   ENDIF
+
+!
+!-- Other Species
+
+                   IF  ( TRIM(emt_att%species_name(ind_inp)) == TRIM(spc_names(ind_mod)) )  THEN
+                      len_index = len_index + 1
+                   ENDIF
+                ENDDO
+             ENDDO
+
+!
+!-- Allocate array for storing the indices of the matched species
+
+             IF  ( len_index > 0 )  THEN
+ 
+                ALLOCATE ( match_spec_input(len_index) ) 
+ 
+                ALLOCATE ( match_spec_model(len_index) )
+
+                IF  ( len_index_voc > 0 )  THEN
+!
+!-- contains indices of the VOC model species
+                   ALLOCATE( match_spec_voc_model(len_index_voc) )
+!
+!-- contains the indices of different values of VOC composition of input variable VOC_composition
+                   ALLOCATE( match_spec_voc_input(len_index_voc) )
+
+                ENDIF
+
+!
+!-- pass the species indices to declared arrays
+
+                len_index = 0
+
+!
+!-- Cycle over model species
+
+                DO ind_mod = 1, nvar
+ 
+!
+!-- Cycle over Input species  
+
+                   DO ind_inp = 1, nspec_emis_inp
+
+!
+!-- VOCs
+
+                      IF  ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND.    &
+                           ALLOCATED(match_spec_voc_input) )  THEN
+                         
+                         DO ind_voc= 1, emt_att%nvoc
+                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                               len_index = len_index + 1
+                               len_index_voc = len_index_voc + 1
+                        
+                               match_spec_input(len_index) = ind_inp
+                               match_spec_model(len_index) = ind_mod
+
+                               match_spec_voc_input(len_index_voc) = ind_voc
+                               match_spec_voc_model(len_index_voc) = ind_mod                         
+                            ENDIF
+                         END DO
+                      ENDIF
+
+!
+!-- Other Species
+
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                         len_index = len_index + 1
+                         match_spec_input(len_index) = ind_inp
+                         match_spec_model(len_index) = ind_mod
+                      ENDIF
+
+                   END DO  ! ind_inp
+                END DO     ! ind_mod
+
+             ELSE  ! if len_index_voc .le. 0
+
+!
+!-- in case there are no species matching
+
+                message_string = 'Non of given emission species'            //         &
+                                 ' matches'                                //          &
+                                 ' model chemical species:'                //          &
+                                 ' Emission routine is not called'  
+                CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
+             ENDIF
+
+!
+!-- Error check (no matching)
+ 
+          ELSE
+
+!
+!-- either spc_names is zero or nspec_emis_inp is not allocated
+             message_string = 'Array of Emission species not allocated:'              // &
+                              ' Either no emission species are provided as input or'  // &
+                              ' no chemical species are used by PALM:'                // &
+                              ' Emission routine is not called'                  
+             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) 
+
+          ENDIF  
+
+!
+!-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
+
+!
+!-- Error check (no species)
+
        CASE DEFAULT
 
-          message_string = 'Emission Module switched ON, but'            //                         &
-                           ' either no emission mode specified or incorrectly given :'  //          &
+          message_string = 'Emission Module switched ON, but'                           //         &
+                           ' either no emission mode specified or incorrectly given :'  //         &
                            ' please, pass the correct value to the namelist parameter "mode_emis"'                  
           CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )              
@@ -701,5 +755,5 @@
        END SELECT
 
-       IF ( debug_output )  CALL debug_message( 'chem_emissions_match', 'end' )
+       IF  ( debug_output )  CALL debug_message( 'chem_emissions_match', 'end' )
 
  END SUBROUTINE chem_emissions_match
@@ -724,5 +778,5 @@
 !    
 !-- Actions for initial runs 
-!  IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+!  IF  (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 !--    ... 
 !    
@@ -735,80 +789,119 @@
 
 
-  IF ( debug_output )  CALL debug_message( 'chem_emissions_init', 'start' )
-
-  !
-  !-- Matching
-  CALL  chem_emissions_match( chem_emis_att, nspec_out )
- 
-  IF ( nspec_out == 0 )  THEN
+    IF  ( debug_output )  CALL debug_message( 'chem_emissions_init', 'start' )
+
+!
+!-- Matching
+
+    CALL  chem_emissions_match( chem_emis_att, n_matched_vars )
+ 
+    IF  ( n_matched_vars == 0 )  THEN
   
-     emission_output_required = .FALSE.
-
-  ELSE
-
-     emission_output_required = .TRUE.
-
-
-     !
-     !-- Set molecule masses' 
-     ALLOCATE( chem_emis_att%xm(nspec_out) )
-
-     DO ispec = 1, nspec_out
-        SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
-           CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2        !< kg/mole
-           CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4        !< kg/mole
-           CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O             !< kg/mole
-           CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2        !< kg/mole   
-           CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3        !< kg/mole
-           CASE ( 'CO'  ); chem_emis_att%xm(ispec) = xm_C + xm_O            !< kg/mole
-           CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2        !< kg/mole
-           CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4        !< kg/mole     
-           CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3  !< kg/mole  
-           CASE DEFAULT
-              chem_emis_att%xm(ispec) = 1.0_wp
-        END SELECT
-     ENDDO
+       emission_output_required = .FALSE.
+
+    ELSE
+
+       emission_output_required = .TRUE.
+
+
+!
+!-- Set molecule masses  (in kg/mol)
+
+       ALLOCATE( chem_emis_att%xm(n_matched_vars) )
+
+       DO ispec = 1, n_matched_vars
+          SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
+             CASE ( 'SO2'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 2
+             CASE ( 'SO4'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 4
+             CASE ( 'NO'   );  chem_emis_att%xm(ispec) = xm_N + xm_O
+             CASE ( 'NO2'  );  chem_emis_att%xm(ispec) = xm_N + xm_O * 2
+             CASE ( 'NH3'  );  chem_emis_att%xm(ispec) = xm_N + xm_H * 3
+             CASE ( 'CO'   );  chem_emis_att%xm(ispec) = xm_C + xm_O
+             CASE ( 'CO2'  );  chem_emis_att%xm(ispec) = xm_C + xm_O * 2
+             CASE ( 'CH4'  );  chem_emis_att%xm(ispec) = xm_C + xm_H * 4
+             CASE ( 'HNO3' );  chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
+             CASE DEFAULT
+                chem_emis_att%xm(ispec) = 1.0_wp
+          END SELECT
+       ENDDO
 
     
-     !
-     !-- assign emission values
-     SELECT CASE ( TRIM( mode_emis ) )   
-
-
-        !
-        !-- PRE-PROCESSED case
-        CASE ( "PRE-PROCESSED" )
-
-           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out) )  
- 
-           !
-           !-- Get emissions at the first time step
-           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
-
-        !
-        !-- Default case
-        CASE ( "DEFAULT" )
-
-           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) ) 
-
-           !
-           !-- Get emissions at the first time step
-           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
-
-        !
-        !-- PARAMETERIZED case
-        CASE ( "PARAMETERIZED" )
-
-           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) )
-
-           !
-           !-- Get emissions at the first time step
-           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out)
-
-     END SELECT
-
-  ENDIF    
-
-  IF ( debug_output )  CALL debug_message( 'chem_emissions_init', 'end' )
+!
+!-- Get emissions for the first time step base on LOD (if defined)
+!-- or emission mode (if no LOD defined)
+
+!
+!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' )
+!--        type of decision structure but I think it is much better
+!--        to implement it this way (i.e., conditional on lod if it
+!--        is defined, and mode if not) as we can easily take out
+!--        the case structure for mode_emis later on.
+!--        (ecc 20140424)
+
+       IF   ( emiss_lod < 0 )  THEN  !-- no LOD defined (not likely)
+
+          SELECT CASE ( TRIM( mode_emis ) )   
+
+             CASE ( 'PARAMETERIZED' )     ! LOD 0
+
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                   ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
+                ENDIF
+
+                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
+
+             CASE ( 'DEFAULT' )           ! LOD 1
+
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                   ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
+                ENDIF
+
+                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
+
+             CASE ( 'PRE-PROCESSED' )     ! LOD 2
+
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                   ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,n_matched_vars) )
+                ENDIF
+ 
+                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
+
+          END SELECT
+
+       ELSE  ! if LOD is defined
+
+          SELECT CASE ( emiss_lod )
+
+             CASE ( 0 )     ! parameterized mode
+
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                   ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
+                ENDIF
+
+                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
+
+             CASE ( 1 )     ! default mode
+
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                   ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
+                ENDIF
+
+                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
+
+             CASE ( 2 )     ! pre-processed mode
+
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                   ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,n_matched_vars) )
+                ENDIF
+ 
+                CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
+
+          END SELECT
+
+       ENDIF
+
+    ENDIF
+
+    IF  ( debug_output )  CALL debug_message( 'chem_emissions_init', 'end' )
 
  END SUBROUTINE chem_emissions_init
@@ -822,5 +915,5 @@
 !-------------------------------------------------------------------------------!
 
- SUBROUTINE chem_emissions_setup( emt_att, emt, nspec_out )
+ SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
  
    USE surface_mod,                                                  &
@@ -835,10 +928,10 @@
   
 
-    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att                         !< variable to store emission information                                                                           
+    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att                         !< variable to store emission information 
 
     TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt  !< variable to store emission input values, 
                                                                                 !< depending on the considered species
 
-    INTEGER,INTENT(IN) ::  nspec_out                                            !< Output of matching routine with number
+    INTEGER,INTENT(IN) ::  n_matched_vars                                       !< Output of matching routine with number
                                                                                 !< of matched species
 
@@ -858,6 +951,6 @@
 
     REAL(wp), DIMENSION(24) :: par_emis_time_factor                             !< time factors for the parameterized mode:
-                                                                               !< fixed houlry profile for example day
-    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  conv_to_ratio          !< factor used for converting input 
+                                                                                !< fixed houlry profile for example day
+    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  conv_to_ratio            !< factor used for converting input 
                                                                                 !< to concentration ratio
     REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  tmp_temp
@@ -883,9 +976,9 @@
     !------------------------------------------------------    
 
-    IF ( emission_output_required )  THEN
+    IF  ( emission_output_required )  THEN
 
        !
        !-- Set emis_dt  
-       IF ( call_chem_at_all_substeps )  THEN
+       IF  ( call_chem_at_all_substeps )  THEN
 
           dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
@@ -897,75 +990,35 @@
        ENDIF
 
-
-       !
-       !-- Conversion of units to the ones employed in PALM 
-       !-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
-
-        IF ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
+ !
+ !-- Conversion of units to the ones employed in PALM 
+ !-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
+
+        IF  ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
 
           SELECT CASE ( TRIM( emt_att%units ) )
-             !
-             !-- kilograms
-             CASE ( 'kg/m2/s', 'KG/M2/S' ) 
-                
-                con_factor=1
-
-             CASE ('kg/m2/hour', 'KG/M2/HOUR' )
-
-                con_factor=hour_to_s
-
-             CASE ( 'kg/m2/day', 'KG/M2/DAY' ) 
-                
-                con_factor=day_to_s
-
-             CASE ( 'kg/m2/year', 'KG/M2/YEAR' )
-              
-                con_factor=year_to_s
-
-             !
-             !-- Tons
-             CASE ( 'ton/m2/s', 'TON/M2/S' ) 
-               
-                con_factor=tons_to_kg
-
-             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ) 
-               
-                con_factor=tons_to_kg*hour_to_s
-
-             CASE ( 'ton/m2/year', 'TON/M2/YEAR' ) 
-               
-                con_factor=tons_to_kg*year_to_s
-
-             !
-             !-- Grams
-             CASE ( 'g/m2/s', 'G/M2/S' )
-          
-                con_factor=g_to_kg
-
-             CASE ( 'g/m2/hour', 'G/M2/HOUR' ) 
-
-                con_factor=g_to_kg*hour_to_s
-
-             CASE ( 'g/m2/year', 'G/M2/YEAR' )
- 
-                con_factor=g_to_kg*year_to_s
-
-             !
-             !-- Micrograms
-             CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' )
- 
-                con_factor=miug_to_kg
-
-             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ) 
- 
-                con_factor=miug_to_kg*hour_to_s
-
-             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' )
-  
-                con_factor=miug_to_kg*year_to_s
+
+             CASE ( 'kg/m2/s',    'KG/M2/S'    );  con_factor = 1.0_wp                   ! kg
+             CASE ( 'kg/m2/hour', 'KG/M2/HOUR' );  con_factor = hour_to_s
+             CASE ( 'kg/m2/day',  'KG/M2/DAY'  );  con_factor = day_to_s
+             CASE ( 'kg/m2/year', 'KG/M2/YEAR' );  con_factor = year_to_s
+
+             CASE ( 'ton/m2/s',    'TON/M2/S'    );  con_factor = tons_to_kg             ! tonnes
+             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' );  con_factor = tons_to_kg*hour_to_s
+             CASE ( 'ton/m2/year', 'TON/M2/YEAR' );  con_factor = tons_to_kg*year_to_s
+
+             CASE ( 'g/m2/s',    'G/M2/S'    );  con_factor = g_to_kg                    ! grams
+             CASE ( 'g/m2/hour', 'G/M2/HOUR' );  con_factor = g_to_kg*hour_to_s
+             CASE ( 'g/m2/year', 'G/M2/YEAR' );  con_factor = g_to_kg*year_to_s
+
+             CASE ( 'micrograms/m2/s',    'MICROGRAMS/M2/S'    );  con_factor = miug_to_kg            ! ug
+             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' );  con_factor = miug_to_kg*hour_to_s
+             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' );  con_factor = miug_to_kg*year_to_s
+
+!
+!-- Error check (need units)
 
              CASE DEFAULT  
-                message_string = 'The Units of the provided emission input' // &
-                                 ' are not the ones required by PALM-4U: please check '      // &
+                message_string = 'The Units of the provided emission input'                 // &
+                                 ' are not the ones required by PALM-4U: please check '     // &
                                  ' emission module documentation.'                                 
                 CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
@@ -973,61 +1026,91 @@
           END SELECT
 
-         
        ENDIF
 
-       !
-       !-- Conversion factor to convert  kg/m**2/s to ppm/s
+!
+!-- Conversion factor to convert  kg/m**2/s to ppm/s
+
        DO i = nxl, nxr
           DO j = nys, nyn
-             !
-             !-- Derive Temperature from Potential Temperature
+
+!
+!-- Derive Temperature from Potential Temperature
+
              tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) * pref_i )**r_cp
-             
-             
-             !>  We need to pass to cssws <- (ppm/s) * dz 
-             !>  Input is Nmole/(m^2*s)
-             !>  To go to ppm*dz multiply the input by (m**2/N)*dz
-             !>  (m**2/N)*dz == V/N 
-             !>  V/N = RT/P 
-             !>  m**3/Nmole               (J/mol)*K^-1           K                      Pa         
+ 
+!            
+!-- We need to pass to cssws <- (ppm/s) * dz 
+!-- Input is Nmole/(m^2*s)
+!-- To go to ppm*dz multiply the input by (m**2/N)*dz
+!-- (m**2/N)*dz == V/N 
+!-- V/N = RT/P 
+
+!                   m**3/Nmole               (J/mol)*K^-1           K                      Pa         
              conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )  
+
           ENDDO
        ENDDO
 
 
-       !
-       !-- Initialize
+!
+!-- Initialize
+
        emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
 
   
-       !
-       !-- PRE-PROCESSED MODE 
-       IF ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
-
-          !
-          !-- Update time indices
+!
+!-- LOD 2 (PRE-PROCESSED MODE)
+
+       IF  ( emiss_lod == 2 )  THEN
+
+! for reference (ecc)
+!       IF  ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
+
+!
+!-- Update time indices
+
           CALL time_preprocessed_indices( index_hh )
 
-       ELSEIF ( TRIM( mode_emis ) == "DEFAULT" )  THEN
-
-          !
-          !-- Allocate array where to store temporary emission values      
-          IF ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
-          !
-          !-- Allocate time factor per category
-          ALLOCATE( time_factor(emt_att%ncat) )  
-          !
-          !-- Read-in hourly emission time factor
-          IF ( TRIM( time_fac_type ) == "HOUR" )  THEN
-
-             !
-             !-- Update time indices
+
+!
+!-- LOD 1 (DEFAULT MODE)
+
+        ELSEIF  ( emiss_lod == 1 )  THEN
+
+! for reference (ecc)
+!       ELSEIF  ( TRIM( mode_emis ) == "DEFAULT" )  THEN
+
+!
+!-- Allocate array where to store temporary emission values
+
+          IF  ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
+
+!
+!-- Allocate time factor per category
+
+          ALLOCATE( time_factor(emt_att%ncat) )
+
+!
+!-- Read-in hourly emission time factor
+
+          IF  ( TRIM(time_fac_type) == "HOUR" )  THEN
+
+!
+!-- Update time indices
+
              CALL time_default_indices( month_of_year, day_of_month, hour_of_day, index_hh )
-             !
-             !-- Check if the index is less or equal to the temporal dimension of HOURLY emission files              
-             IF ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) )  THEN
-                !
-                !-- Read-in the correspondant time factor              
+
+!
+!-- Check if the index is less or equal to the temporal dimension of HOURLY emission files
+
+             IF  ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) )  THEN
+
+!
+!-- Read-in the correspondant time factor
+
                 time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh)      
+
+!
+!-- Error check (time out of range)
 
              ELSE
@@ -1038,20 +1121,29 @@
 
              ENDIF
-          !
-          !-- Read-in MDH emissions time factors
-          ELSEIF ( TRIM( time_fac_type ) == "MDH" )  THEN
-
-             !
-             !-- Update time indices      
+
+!
+!-- Read-in MDH emissions time factors
+
+          ELSEIF  ( TRIM( time_fac_type ) == "MDH" )  THEN
+
+!
+!-- Update time indices
              CALL time_default_indices( daytype_mdh, month_of_year, day_of_month,     &
                   hour_of_day, index_mm, index_dd,index_hh )
 
-             !
-             !-- Check if the index is less or equal to the temporal dimension of MDH emission files              
-             IF ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )  THEN
-                !
-                !-- Read-in the correspondant time factor              
-                time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) * emt_att%mdh_emis_time_factor(:,index_dd) *   &
-                                                                         emt_att%mdh_emis_time_factor(:,index_hh)
+!
+!-- Check if the index is less or equal to the temporal dimension of MDH emission files
+
+             IF  ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )  THEN
+
+!
+!-- Read in corresponding time factor
+
+                time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) *    &
+                                 emt_att%mdh_emis_time_factor(:,index_dd) *    &
+                                 emt_att%mdh_emis_time_factor(:,index_hh)
+
+!
+!-- Error check (MDH time factor not provided)
      
              ELSE
@@ -1061,5 +1153,8 @@
                 CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
 
-             ENDIF 
+             ENDIF  
+
+!
+!-- Error check (no time factor defined)
 
           ELSE
@@ -1071,177 +1166,212 @@
           ENDIF
 
-       !
-       !-- PARAMETERIZED MODE
-       ELSEIF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
+!
+!-- PARAMETERIZED MODE
+
+       ELSEIF ( emiss_lod == 0 )  THEN
+
+
+! for reference (ecc)
+!       ELSEIF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
         
-          !
-          !-- assign constant values of time factors, diurnal time profile for traffic sector
-          par_emis_time_factor( : ) =  &
-            (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, 0.056, 0.053, 0.051, 0.051, 0.052, 0.055,  &
-               0.059, 0.061, 0.064, 0.067, 0.069, 0.069, 0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
+!
+!-- assign constant values of time factors, diurnal profile for traffic sector
+
+          par_emis_time_factor( : ) =  (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039,      &
+                                          0.056, 0.053, 0.051, 0.051, 0.052, 0.055,      &
+                                          0.059, 0.061, 0.064, 0.067, 0.069, 0.069,      &
+                                          0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
           
-          IF ( .NOT. ALLOCATED( time_factor ) ) ALLOCATE( time_factor(1) )
-
-          !
-          !-- Get time-factor for specific hour
+          IF  ( .NOT. ALLOCATED (time_factor) )  ALLOCATE (time_factor(1))
+
+!
+!-- Get time-factor for specific hour
+
           index_hh = hour_of_day
-
           time_factor(1) = par_emis_time_factor(index_hh)
 
-       ENDIF
+       ENDIF  ! emiss_lod
 
        
-       !
-       !--  Emission distribution calculation
-
-       !
-       !-- PARAMETERIZED case
-       IF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
-
-          DO ispec = 1, nspec_out
-
-             !
-             !-- Units are micromoles/m**2*day (or kilograms/m**2*day for PMs)
-             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = surface_csflux(match_spec_input(ispec)) *  &
-                                                           time_factor(1) * hour_to_s
+!
+!--  Emission distribution calculation
+
+!
+!-- LOD 0 (PARAMETERIZED mode)
+
+       IF  ( emiss_lod == 0 )  THEN
+
+! for reference (ecc)
+!       IF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
+
+          DO ispec = 1, n_matched_vars
+
+!
+!-- Units are micromoles/m**2*day (or kilograms/m**2*day for PMs)
+
+             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =            &
+                       surface_csflux(match_spec_input(ispec)) *     &
+                       time_factor(1) * hour_to_s
 
           ENDDO 
 
-       !
-       !-- PRE-PROCESSED case
-       ELSEIF ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
-
-          !
-          !-- Cycle over species:
-          !-- nspec_out represents the number of species in common between the emission input data
-          !-- and the chemistry mechanism used
-          DO ispec=1,nspec_out  
-    
-             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%                                 &
-                                                            preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1) * &
-                                                                con_factor
           
-          ENDDO
-          
-       !
-       !-- DEFAULT case
-       ELSEIF ( TRIM( mode_emis ) == "DEFAULT" )  THEN
-
-          !
-          !-- Allocate array for the emission value corresponding to a specific category and time factor
-          ALLOCATE( delta_emis(nys:nyn,nxl:nxr) )  
-
-          !
-          !-- Cycle over categories
+!
+!-- LOD 1 (DEFAULT mode)
+
+
+       ELSEIF  ( emiss_lod == 1 )  THEN
+
+! for referene (ecc)
+!       ELSEIF  ( TRIM( mode_emis ) == "DEFAULT" )  THEN
+
+!
+!-- Allocate array for the emission value corresponding to a specific category and time factor
+
+          ALLOCATE (delta_emis(nys:nyn,nxl:nxr))
+
+!
+!-- Cycle over categories
+
           DO icat = 1, emt_att%ncat
  
-             !
-             !-- Cycle over Species
-             !-- nspec_out represents the number of species in common between the emission input data
-             !-- and the chemistry mechanism used
-             DO ispec = 1, nspec_out
-
-                emis(nys:nyn,nxl:nxr) = emt(match_spec_input(ispec))%default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
-
-
-                !
-                !-- NOx
-                IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "NO" )  THEN
+!
+!-- Cycle over Species:  n_matched_vars represents the number of species
+!-- in common between the emission input data and the chemistry mechanism used
+
+             DO ispec = 1, n_matched_vars
+
+                emis(nys:nyn,nxl:nxr) =                    &
+                       emt(match_spec_input(ispec))%       &
+                           default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
+
+!
+!-- NO 
+
+                IF  ( TRIM(spc_names(match_spec_model(ispec))) == "NO" )  THEN
                 
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
-                                                  emt_att%nox_comp(icat,1) * con_factor * hour_per_day
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                 delta_emis(nys:nyn,nxl:nxr)
-
-                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "NO2" )  THEN
-
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
-                                                  emt_att%nox_comp(icat,2) * con_factor * hour_per_day
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                 delta_emis(nys:nyn,nxl:nxr)
- 
-                !
-                !-- SOx
-                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "SO2" )  THEN
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
+                                                 time_factor(icat) *           &
+                                                 emt_att%nox_comp(icat,1) *    &
+                                                 con_factor * hour_per_day
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
+                               delta_emis(nys:nyn,nxl:nxr)
+!
+!-- NO2
+
+                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "NO2" )  THEN
+
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
+                                                 time_factor(icat) *           &
+                                                 emt_att%nox_comp(icat,2) *    &
+                                                 con_factor * hour_per_day
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
+                               delta_emis(nys:nyn,nxl:nxr)
+ 
+!
+!-- SO2
+                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "SO2" )  THEN
                   
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
-                                                 emt_att%sox_comp(icat,1) * con_factor * hour_per_day
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                 delta_emis(nys:nyn,nxl:nxr)
-
-                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "SO4" )  THEN
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
+                                                 time_factor(icat) *           &
+                                                 emt_att%sox_comp(icat,1) *    &
+                                                 con_factor * hour_per_day
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
+                               delta_emis(nys:nyn,nxl:nxr)
+
+!
+!-- SO4
+
+                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "SO4" )  THEN
+
                   
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &                 !<  kg/m2*s
-                                                 emt_att%sox_comp(icat,2) * con_factor * hour_per_day
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                 delta_emis(nys:nyn,nxl:nxr)
- 
-
-                !
-                !-- PMs
-                !-- PM1
-                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1" )  THEN
-
-                   !
-                   !-- Cycle over PM1 components
-                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,1) )
-
-                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
-                                                     emt_att%pm_comp(icat,i_pm_comp,1) * con_factor * hour_per_day 
-
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                    delta_emis(nys:nyn,nxl:nxr)
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
+                                                 time_factor(icat) *           &
+                                                 emt_att%sox_comp(icat,2) *    &
+                                                 con_factor * hour_per_day
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
+                               delta_emis(nys:nyn,nxl:nxr)
+ 
+
+!
+!-- PM1
+
+                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1" )  THEN
+
+                   DO  i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,1) )   ! cycle through components
+
+                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *              &     ! kg/m2/s
+                                                    time_factor(icat) *                  &
+                                                    emt_att%pm_comp(icat,i_pm_comp,1) *  &
+                                                    con_factor * hour_per_day 
+
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
+                               delta_emis(nys:nyn,nxl:nxr)
+
                    ENDDO
 
-                !
-                !-- PM2.5
-                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM25" )  THEN
-
-                   !
-                   !-- Cycle over PM2.5 components
-                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,2) )
-
-                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
-                                                     emt_att%pm_comp(icat,i_pm_comp,2) * con_factor * hour_per_day 
-
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                    delta_emis(nys:nyn,nxl:nxr)
+!
+!-- PM2.5
+
+                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM25" )  THEN
+
+                   DO  i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,2) )   ! cycle through components
+
+                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *              &     ! kg/m2/s
+                                                    time_factor(icat) *                  &
+                                                    emt_att%pm_comp(icat,i_pm_comp,2) *  &
+                                                    con_factor * hour_per_day 
+
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
+                               delta_emis(nys:nyn,nxl:nxr)
  
                    ENDDO
 
-                !
-                !-- PM10 
-                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
-
-                   !
-                   !-- Cycle over PM10 components
-                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) )  
-
-                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
-                                                     emt_att%pm_comp(icat,i_pm_comp,3) * con_factor * hour_per_day 
-
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+ &
-                                                                 delta_emis(nys:nyn,nxl:nxr) 
+!
+!-- PM10
+
+                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM10" )  THEN
+
+                   DO  i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) )   ! cycle through components
+
+                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *              &     ! kg/m2/s
+                                                    time_factor(icat)     *              &
+                                                    emt_att%pm_comp(icat,i_pm_comp,3) *  &
+                                                    con_factor * hour_per_day 
+
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
+                               delta_emis(nys:nyn,nxl:nxr) 
 
                    ENDDO
 
-                !
-                !-- VOCs
-                ELSEIF ( SIZE( match_spec_voc_input ) > 0 )  THEN
-
-                   DO ivoc = 1, SIZE( match_spec_voc_input )
-
-                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == TRIM( emt_att%voc_name(ivoc) ) )  THEN    
-
-                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
-                                                        emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) * &
-                                                         con_factor * hour_per_day
-
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                       delta_emis(nys:nyn,nxl:nxr) 
+!
+!-- VOCs
+
+                ELSEIF  ( SIZE( match_spec_voc_input ) > 0 )  THEN
+
+                   DO  ivoc = 1, SIZE( match_spec_voc_input )          ! cycle through components
+
+                      IF  ( TRIM(spc_names(match_spec_model(ispec))) ==                &
+                            TRIM(emt_att%voc_name(ivoc)) )  THEN    
+
+                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *           &
+                                                       time_factor(icat) *               &
+                                                       emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) *   &
+                                                       con_factor * hour_per_day
+
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                    &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
+                               delta_emis(nys:nyn,nxl:nxr) 
 
                       ENDIF                       
@@ -1249,15 +1379,18 @@
                    ENDDO
                 
-                !
-                !-- any other species
+!
+!-- any other species
+
                 ELSE
 
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
-                                                  con_factor * hour_per_day
- 
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
-                                                                 delta_emis(nys:nyn,nxl:nxr) 
-
-                ENDIF
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *                 &
+                                                 time_factor(icat) *                     &
+                                                 con_factor * hour_per_day
+ 
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                          &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +              &
+                               delta_emis(nys:nyn,nxl:nxr) 
+
+                ENDIF  ! TRIM spc_names
                 
                 emis(:,:)= 0 
@@ -1269,120 +1402,173 @@
           ENDDO
 
-       ENDIF
+!
+!-- LOD 2 (PRE-PROCESSED mode)
+
+       ELSEIF  ( emiss_lod == 2 )  THEN
+
+! for reference (ecc)
+!       ELSEIF  ( TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
+
+!
+!-- Cycle over species: n_matched_vars represents the number of species
+!-- in common between the emission input data and the chemistry mechanism used
+
+          DO ispec = 1, n_matched_vars  
+ 
+! (ecc)   
+             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                &
+                       emt(match_spec_input(ispec))%                                     &
+                           preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1) *   &
+                       con_factor
+
+
+!             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                &
+!                       emt(match_spec_input(ispec))%                                     &
+!                           preproc_emission_data(index_hh,1,:,:) *   &
+!                       con_factor
+          ENDDO
+
+       ENDIF  ! emiss_lod
 
        
 !
 !-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
-!
-!-- PARAMETERIZED mode
-       !
-       !-- Units of inputs are micromoles/(m**2*s)
-       IF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
-
-          IF ( street_type_f%from_file )  THEN
-
-!
-!--          Streets are lsm surfaces, hence, no usm surface treatment required. 
-!--          However, urban surface may be initialized via default initialization
-!--          in surface_mod, e.g. at horizontal urban walls that are at k == 0
-!--          (building is lower than the first grid point). Hence, in order to 
-!--          have only emissions at streets, set the surfaces emissions to zero
-!--          at urban walls. 
-             IF ( surf_usm_h%ns >=1 )  surf_usm_h%cssws = 0.0_wp
-!
-!--          Now, treat land-surfaces.
+
+!
+!-- LOD 0 (PARAMETERIZED mode)
+!-- Units of inputs are micromoles/m2/s
+
+       IF  ( emiss_lod == 0 )  THEN
+! for reference (ecc)
+!       IF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
+
+          IF  (street_type_f%from_file)  THEN
+
+!
+!-- Streets are lsm surfaces, hence, no usm surface treatment required. 
+!-- However, urban surface may be initialized via default initialization
+!-- in surface_mod, e.g. at horizontal urban walls that are at k == 0
+!-- (building is lower than the first grid point). Hence, in order to 
+!-- have only emissions at streets, set the surfaces emissions to zero
+!-- at urban walls.
+ 
+             IF  ( surf_usm_h%ns >=1 )  surf_usm_h%cssws = 0.0_wp
+
+!
+!-- Treat land-surfaces.
+
              DO  m = 1, surf_lsm_h%ns
+
                 i = surf_lsm_h%i(m)
                 j = surf_lsm_h%j(m)
                 k = surf_lsm_h%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 matched species
-                   DO  ispec = 1, nspec_out
-
-                      !
-                      !-- PMs are already in 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
-
-                         !
-                         !--           kg/(m^2*s) * kg/m^3
-                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
-                                                                        emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
-                                                                         rho_air(k)                                    !< in kg/m^3
+                IF  ( street_type_f%var(j,i) >= main_street_id    .AND.                  &
+                      street_type_f%var(j,i) < max_street_id   )  THEN
+
+!
+!-- Cycle over matched species
+
+                   DO  ispec = 1, n_matched_vars
+
+!
+!-- PMs are already in 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
+
+!
+!-- kg/(m^2*s) * kg/m^3
+                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
+                                   emiss_factor_main(match_spec_input(ispec)) *          &
+                                   emis_distribution(1,j,i,ispec) *                      &     ! kg/(m^2*s)
+                                   rho_air(k)                                                  ! kg/m^3
                          
-                      !
-                      !-- Other Species
-                      !-- Inputs are micromoles
+!
+!-- Other Species
+!-- Inputs are micromoles
+
                       ELSE
 
-                         !   
-                         !--             ppm/s *m *kg/m^3                
-                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
-                                                                        emis_distribution(1,j,i,ispec) *            &  !< in micromoles/(m^2*s)
-                                                                         conv_to_ratio(k,j,i) *                     &  !< in m^3/Nmole     
-                                                                          rho_air(k)                                   !< in kg/m^3
+!   
+!-- ppm/s *m *kg/m^3                
+                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
+                                   emiss_factor_main(match_spec_input(ispec)) *          &
+                                   emis_distribution(1,j,i,ispec) *                      &     ! micromoles/(m^2*s)
+                                   conv_to_ratio(k,j,i) *                                &     ! m^3/Nmole     
+                                   rho_air(k)                                                  ! kg/m^3
 
                       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 matched species
-                   DO  ispec = 1, nspec_out
-
-                      !
-                      !-- PMs are already in 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
-
-                         !
-                         !--           kg/(m^2*s) * kg/m^3
-                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
-                                                                        emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
-                                                                         rho_air(k)                                    !< in kg/m^3   
-                      !
-                      !-- Other species
-                      !-- Inputs are micromoles
+                   ENDDO  ! ispec
+
+
+                ELSEIF  ( street_type_f%var(j,i) >= side_street_id   .AND.               &
+                          street_type_f%var(j,i) < main_street_id )  THEN
+
+!
+!-- Cycle over matched species
+
+                   DO  ispec = 1, n_matched_vars
+
+!
+!-- PMs are already in 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
+
+!
+!-- kg/(m^2*s) * kg/m^3
+                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
+                                   emiss_factor_side(match_spec_input(ispec)) *          &
+                                   emis_distribution(1,j,i,ispec) *                      &     ! kg/(m^2*s)
+                                   rho_air(k)                                                  ! kg/m^3   
+!
+!-- Other species
+!-- Inputs are micromoles
+
                       ELSE
-                    
-                         !   
-                         !--             ppm/s *m *kg/m^3      
-                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
-                                                                        emis_distribution(1,j,i,ispec) *            &  !< in micromoles/(m^2*s)
-                                                                         conv_to_ratio(k,j,i) *                     &  !< in m^3/Nmole       
-                                                                          rho_air(k)                                   !< in kg/m^3   
+
+!   
+!-- ppm/s *m *kg/m^3
+
+                         surf_lsm_h%cssws(match_spec_model(ispec),m) =                   &
+                                   emiss_factor_side(match_spec_input(ispec)) *          &
+                                   emis_distribution(1,j,i,ispec) *                      &  ! micromoles/(m^2*s)
+                                   conv_to_ratio(k,j,i) *                                &  ! m^3/Nmole       
+                                   rho_air(k)                                               ! kg/m^3   
+
                       ENDIF
 
-                   ENDDO
+                   ENDDO  ! ispec
+
+!
+!-- If no street type is defined, then assign zero emission to all the species
 
                 ELSE
 
-                   !
-                   !-- If no street type is defined, then assign zero emission to all the species
                    surf_lsm_h%cssws(:,m) = 0.0_wp
 
-                ENDIF
-
-             ENDDO
-
-          ENDIF
-
-       !
-       !-- For both DEFAULT and PRE-PROCESSED mode
+                ENDIF  ! street type
+
+             ENDDO  ! m
+
+          ENDIF  ! street_type_f%from_file
+
+
+!
+!-- LOD 1 (DEFAULT) and LOD 2 (PRE-PROCESSED)
+
+
        ELSE    
        
 
-          DO ispec = 1, nspec_out 
+          DO ispec = 1, n_matched_vars 
                    
 !
 !--          Default surfaces
+
              DO  m = 1, surf_def_h(0)%ns
 
@@ -1390,55 +1576,54 @@
                 j = surf_def_h(0)%j(m)
 
-                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
-
-                   !
-                   !-- PMs
-                   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
+                IF  ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
+
+!
+!-- PMs
+                   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
                 
-                      !
-                      !--            kg/(m^2*s) *kg/m^3                         kg/(m^2*s)                 
-                      surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*  &
-                                                                        rho_air(nzb)                           !< in kg/m^3  
- 
+                      surf_def_h(0)%cssws(match_spec_model(ispec),m) =         &     ! kg/m2/s * kg/m3
+                               emis_distribution(1,j,i,ispec)*                 &     ! kg/m2/s
+                               rho_air(nzb)                                          ! kg/m^3  
  
                    ELSE
 
-                      !
-                      !-- VOCs
-                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
-                         !
-                         !--           (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
-                         surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
-                                                                           conv_to_ratio(nzb,j,i) *         &  !< in m^3/mole
-                                                                            ratio2ppm *                     &  !< in ppm
-                                                                             rho_air(nzb)                      !< in kg/m^3  
-
-
-                      !
-                      !-- Other species
+!
+!-- VOCs
+                      IF  ( len_index_voc > 0                                           .AND.      &
+                            emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
+
+                         surf_def_h(0)%cssws(match_spec_model(ispec),m) =      &  ! ppm/s * m * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &  ! mole/m2/s
+                               conv_to_ratio(nzb,j,i) *                        &  ! m^3/mole
+                               ratio2ppm *                                     &  ! ppm
+                               rho_air(nzb)                                       ! kg/m^3  
+
+
+!
+!-- Other species
+
                       ELSE
 
-                         !
-                         !--           (ppm/s) * m  * kg/m^3                    kg/(m^2/s)
-                         surf_def_h(0)%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &       
-                                                                           ( 1.0_wp / emt_att%xm(ispec) ) * &  !< in mole/kg 
-                                                                             conv_to_ratio(nzb,j,i) *       &  !< in m^3/mole  
-                                                                              ratio2ppm *                   &  !< in ppm 
-                                                                               rho_air(nzb)                    !< in  kg/m^3  
+                         surf_def_h(0)%cssws(match_spec_model(ispec),m) =      &   ! ppm/s * m * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &   ! kg/m2/s
+                               ( 1.0_wp / emt_att%xm(ispec) ) *                &   ! mole/kg 
+                               conv_to_ratio(nzb,j,i) *                        &   ! m^3/mole  
+                               ratio2ppm *                                     &   ! ppm 
+                               rho_air(nzb)                                        !  kg/m^3  
   
-
-                      ENDIF
-
-                   ENDIF
-
-                ENDIF
-
-             ENDDO
+                      ENDIF  ! VOC
+
+                   ENDIF  ! PM
+
+                ENDIF  ! emis_distribution > 0
+
+             ENDDO  ! m
          
 !
-!--          LSM surfaces
-             DO m = 1, surf_lsm_h%ns
+!-- LSM surfaces
+
+             DO  m = 1, surf_lsm_h%ns
 
                 i = surf_lsm_h%i(m)
@@ -1446,51 +1631,54 @@
                 k = surf_lsm_h%k(m)
 
-                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
-
-                   !
-                   !-- PMs
-                   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
-   
-                      !
-                      !--         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)            
-                      surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
-                                                                     rho_air(k)                                !< in kg/m^3
+                IF  ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
+
+!
+!-- PMs
+                   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_lsm_h%cssws(match_spec_model(ispec),m) =            &    ! kg/m2/s * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &    ! kg/m2/s
+                               rho_air(k)                                           ! kg/m^3
  
                    ELSE
 
-                      !
-                      !-- VOCs
-                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
-                         !
-                         !--          (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
-                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &     
-                                                                        conv_to_ratio(k,j,i) *           &     !< in m^3/mole
-                                                                         ratio2ppm *                     &     !< in ppm
-                                                                          rho_air(k)                           !< in kg/m^3  
-
-                      !
-                      !-- Other species
+!
+!-- VOCs
+
+                      IF  ( len_index_voc > 0                                           .AND.      &
+                            emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
+
+                         surf_lsm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &   ! mole/m2/s  
+                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
+                               ratio2ppm *                                     &   ! ppm
+                               rho_air(k)                                          ! kg/m^3  
+
+!
+!-- Other species
+
                       ELSE
-                         !
-                         !--           (ppm/s) * m  * kg/m^3                    kg/(m^2/s)
-                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  & 
-                                                                        ( 1.0_wp / emt_att%xm(ispec) ) * &     !< in mole/kg
-                                                                         conv_to_ratio(k,j,i) *          &     !< in m^3/mole
-                                                                          ratio2ppm *                    &     !< in ppm
-                                                                           rho_air(k)                          !< in kg/m^3     
+
+                         surf_lsm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &   ! kg/m2/s
+                               ( 1.0_wp / emt_att%xm(ispec) ) *                &   ! mole/kg
+                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
+                               ratio2ppm *                                     &   ! ppm
+                               rho_air(k)                                          ! kg/m^3     
                                                    
-                      ENDIF
-
-                   ENDIF
-
-                ENDIF
-
-             ENDDO
-
-!
-!--          USM surfaces
-             DO m = 1, surf_usm_h%ns
+                      ENDIF  ! VOC
+
+                   ENDIF  ! PM
+
+                ENDIF  ! emis_distribution
+
+             ENDDO  ! m
+
+!
+!-- USM surfaces
+
+             DO  m = 1, surf_usm_h%ns
 
                 i = surf_usm_h%i(m)
@@ -1498,50 +1686,48 @@
                 k = surf_usm_h%k(m)
 
-                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
-
-                   !
-                   !-- PMs
-                   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
+                IF  ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
+
+!
+!-- PMs
+                   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
                    
-                      !
-                      !--          kg/(m^2*s) *kg/m^3                             kg/(m^2*s)                     
-                      surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec)*  & 
-                                                                     rho_air(k)                                !< in kg/m^3
-
- 
+                      surf_usm_h%cssws(match_spec_model(ispec),m) =            &    ! kg/m2/s * kg/m3
+                               emis_distribution(1,j,i,ispec)*                 &    ! kg/m2/s
+                               rho_air(k)                                           ! kg/m^3
+
                    ELSE
                       
-                      !
-                      !-- VOCs
-                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
-                         !
-                         !--          (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
-                         surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
-                                                                        conv_to_ratio(k,j,i) *           &     !< in m^3/mole
-                                                                         ratio2ppm *                     &     !< in ppm
-                                                                          rho_air(k)                           !< in kg/m^3   
-
-                      !
-                      !-- Other species
+!
+!-- VOCs
+                      IF  ( len_index_voc > 0                                           .AND.      &
+                            emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
+
+                         surf_usm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &   ! m2/s
+                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
+                               ratio2ppm *                                     &   ! ppm
+                               rho_air(k)                                          ! kg/m^3   
+
+!
+!-- Other species
                       ELSE
 
-                         !
-                         !--          (ppm/s) * m * kg/m^3                        kg/(m^2/s)                      
-                         surf_usm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &   
-                                                                        ( 1.0_wp / emt_att%xm(ispec) ) * &     !< in mole/kg
-                                                                         conv_to_ratio(k,j,i) *          &     !< in m^3/mole
-                                                                          ratio2ppm*                     &     !< in ppm
-                                                                           rho_air(k)                          !< in kg/m^3   
-
-
-                      ENDIF
-
-                   ENDIF
-
-                ENDIF
- 
-             ENDDO
+                         surf_usm_h%cssws(match_spec_model(ispec),m) =         &   ! ppm/s * m * kg/m3
+                               emis_distribution(1,j,i,ispec) *                &   ! kg/m2/s
+                               ( 1.0_wp / emt_att%xm(ispec) ) *                &   ! mole/kg
+                               conv_to_ratio(k,j,i) *                          &   ! m^3/mole
+                               ratio2ppm*                                      &   ! ppm
+                               rho_air(k)                                          ! kg/m^3   
+
+
+                      ENDIF  ! VOC
+
+                   ENDIF  ! PM
+
+                ENDIF  ! emis_distribution
+ 
+             ENDDO  ! m
 
           ENDDO
@@ -1549,7 +1735,8 @@
        ENDIF
 
-       !
-       !-- Deallocate time_factor in case of DEFAULT mode) 
-       IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor )
+!
+!-- Deallocate time_factor in case of DEFAULT mode) 
+
+       IF  ( ALLOCATED (time_factor) )  DEALLOCATE (time_factor)
 
    ENDIF
