Index: /palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3967)
+++ /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
Index: /palm/trunk/SOURCE/chem_modules.f90
===================================================================
--- /palm/trunk/SOURCE/chem_modules.f90	(revision 3967)
+++ /palm/trunk/SOURCE/chem_modules.f90	(revision 3968)
@@ -27,4 +27,13 @@
 ! -----------------
 ! $Id$
+! - introduced namelist item chem_modules@emiss_lod as future
+! - replacement to chem_modules@mode_emis.  Currently keeping both
+!   for backward compatibility.  chem_modules@mode_emis will be
+!   depreciated upon migration of all dependent modules (e.g., salsa)
+!   to chem_modules@emiss_lod
+!
+! (ecc) 20190513 replaced nspec_out with n_matched_vars
+! 
+! 3877 2019-04-08 19:09:16Z knoop
 ! Formatting, clean-up, clarified/corrected comments
 ! 
@@ -102,4 +111,6 @@
     INTEGER(iwp) ::  cs_pr_count                           = 0      !< counter for chemical species profiles
     INTEGER(iwp) ::  cs_vertical_gradient_level_ind(99,10) = -9999  !< grid index values of cs_vertical_gradient_level
+    INTEGER(iwp) ::  emiss_lod                             = -1     !< namelist parameter: chem emission LOD (same as mode_emis)
+                                                                    !< -1 = unassigned, 0 = parameterized, 1 = default, 2 = pre-processed
     INTEGER(iwp) ::  ibc_cs_b                                       !< integer flag for bc_cs_b
     INTEGER(iwp) ::  ibc_cs_t                                       !< integer flag for bc_cs_t
@@ -107,5 +118,5 @@
     INTEGER(iwp) ::  max_pr_cs                             = 0      !< 
     INTEGER(iwp) ::  max_street_id                         = 0      !< namelist parameter: upper bound of main street IDs (OpenStreetMaps) for PARAMETERIZED mode      
-    INTEGER(iwp) ::  nspec_out                                      !< output of routine chem_emis_matching with.....???
+    INTEGER(iwp) ::  n_matched_vars                                 !< number of matched emissions variables
     INTEGER(iwp) ::  side_street_id                        = 0      !< namelist parameter: lower bound of side street IDs (OpenStreetMaps) for PARAMETERIZED mode
 
Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3967)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3968)
@@ -27,4 +27,13 @@
 ! -----------------
 ! $Id$
+! - added "emiss_lod" which serves the same function as "mode_emis"
+!   both will be synchronized with emiss_lod having pirority over
+!   mode_emis (see informational messages)
+! - modified existing error message and introduced new informational messages
+!    - CM0436 - now also applies to invalid LOD settings
+!    - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
+!    - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
+! 
+! 3930 2019-04-24 14:57:18Z forkel
 ! Changed chem_non_transport_physics to chem_non_advective_processes
 ! 
@@ -1682,11 +1691,14 @@
 !
 !-- Emission mode info
-    IF ( mode_emis == "DEFAULT" )  THEN 
-       WRITE( io, 5 ) 
-    ELSEIF ( mode_emis == "PARAMETERIZED" )  THEN 
-       WRITE( io, 6 )
-    ELSEIF ( mode_emis == "PRE-PROCESSED" )  THEN 
-       WRITE( io, 7 )
-    ENDIF 
+!-- At the moment the evaluation is done with both emiss_lod and mode_emis
+!-- but once salsa has been migrated to emiss_lod the .OR. mode_emis
+!-- conditions can be removed (ecc 20190513)
+    IF     ( (emiss_lod == 1) .OR. (mode_emis == 'DEFAULT') )        THEN
+        WRITE ( io, 5 )
+    ELSEIF ( (emiss_lod == 0) .OR. (mode_emis == 'PARAMETERIZED') )  THEN 
+        WRITE ( io, 6 )
+    ELSEIF ( (emiss_lod == 2) .OR. (mode_emis == 'PRE-PROCESSED') )  THEN 
+        WRITE ( io, 7 )
+    ENDIF
 !
 !-- Photolysis scheme info
@@ -2209,5 +2221,6 @@
          decycle_method,                   &
          deposition_dry,                   &
-         emissions_anthropogenic,          &  
+         emissions_anthropogenic,          &
+         emiss_lod,                        &
          emiss_factor_main,                &
          emiss_factor_side,                &                      
@@ -2271,7 +2284,67 @@
 !
 !--    check for emission mode for chem species
-    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
-       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
-       CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
+
+    IF ( emiss_lod < 0 )  THEN   !- if LOD not defined in namelist
+       IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.      &
+            ( mode_emis /= 'DEFAULT'        )    .AND.      &
+            ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
+          message_string = 'Incorrect mode_emiss  option select. Please check spelling'
+          CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
+       ENDIF
+    ELSE
+       IF ( ( emiss_lod /= 0 )    .AND.         &
+            ( emiss_lod /= 1 )    .AND.         &
+            ( emiss_lod /= 2 ) )  THEN
+          message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
+          CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
+       ENDIF
+    ENDIF
+
+! for reference (ecc)
+!    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
+!       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
+!       CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
+!    ENDIF
+
+!
+!-- conflict resolution for emiss_lod and mode_emis
+!-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
+!-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
+!-- this check is in place to retain backward compatibility with salsa until the
+!-- code is migrated completed to emiss_lod
+    IF  ( emiss_lod >= 0 ) THEN
+       SELECT CASE  ( emiss_lod )
+          CASE (0)  !- parameterized mode
+             mode_emis = 'PARAMETERIZED'
+          CASE (1)  !- default mode
+             mode_emis = 'DEFAULT'
+          CASE (2)  !- preprocessed mode
+             mode_emis = 'PRE-PROCESSED'
+       END SELECT
+       
+       message_string = 'Synchronizing mode_emis to defined emiss_lod'               //  &
+                        CHAR(10)  //  '          '                                   //  &
+                        'NOTE - mode_emis will be depreciated in future releases'    //  &
+                        CHAR(10)  //  '          '                                   //  &
+                        'please use emiss_lod to define emission mode'
+
+       CALL message ( 'parin_chem', 'CM0463', 0, 0, 0, 6, 0 )
+    ELSE  ! if emiss_lod is not set
+       SELECT CASE ( mode_emis )
+          CASE ('PARAMETERIZED')
+             emiss_lod = 0
+          CASE ('DEFAULT')
+             emiss_lod = 1
+          CASE ('PRE-PROCESSED')
+             emiss_lod = 2
+       END SELECT
+
+       message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //  &
+                        CHAR(10)  //  '          '                                   //  &
+                        'NOTE - mode_emis will be depreciated in future releases'    //  &
+                        CHAR(10)  //  '          '                                   //  &
+                        '       please use emiss_lod to define emission mode'
+
+       CALL message ( 'parin_chem', 'CM0464', 0, 0, 0, 6, 0 )
     ENDIF
 
Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3967)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3968)
@@ -25,4 +25,18 @@
 ! -----------------
 ! $Id$
+! - clean-up index notations for emission_values to eliminate magic numbers
+! - introduce temporary variable dum_var_5d as well as subroutines
+!   get_var_5d_real and get_var_5d_real_dynamic
+! - remove emission-specific code in generic get_variable routines
+! - in subroutine netcdf_data_input_chemistry_data change netCDF LOD 1 
+!   (default) emission_values to the following index order:
+!   z, y, x, species, category
+! - in subroutine netcdf_data_input_chemistry_data
+!   changed netCDF LOD 2 pre-processed emission_values to the following index 
+!   order: time, z, y, x, species
+! - in type chem_emis_att_type replace nspec with n_emiss_species
+!   but retained nspec for backward compatibility with salsa_mod. (E.C. Chan)
+! 
+! 3961 2019-05-08 16:12:31Z suehring
 ! Revise checks for building IDs and types
 ! 
@@ -289,10 +303,9 @@
 ! Initial revision (suehring)
 !
-!
-!
-!
 ! Authors:
 ! --------
 ! @author Matthias Suehring
+! @author Edward C. Chan
+! @author Emanuele Russo
 !
 ! Description:
@@ -300,10 +313,17 @@
 !> Modulue contains routines to input data according to Palm input data
 !> standart using dynamic and static input files.
-!> @todo - Chemistry: revise reading of netcdf file and ajdust formatting according to standard!!!
+!> @todo - Chemistry: revise reading of netcdf file and ajdust formatting
+!>         according to standard!!! (ecc/done)
 !> @todo - Order input alphabetically
 !> @todo - Revise error messages and error numbers
 !> @todo - Input of missing quantities (chemical species, emission rates)
 !> @todo - Defninition and input of still missing variable attributes
+!>         (ecc/what are they?)
 !> @todo - Input of initial geostrophic wind profiles with cyclic conditions.
+!> @todo - remove z dimension from default_emission_data nad preproc_emission_data
+!          and correpsonding subroutines get_var_5d_real and get_var_5d_dynamic (ecc)
+!> @todo - decpreciate chem_emis_att_type@nspec (ecc)
+!> @todo - depreciate subroutines get_variable_4d_to_3d_real and
+!>         get_variable_5d_to_4d_real (ecc)
 !------------------------------------------------------------------------------!
  MODULE netcdf_data_input_mod
@@ -348,5 +368,5 @@
        CHARACTER(LEN=17) ::  char_s = 'ls_forcing_south_' !< leading substring for variables at south boundary 
        CHARACTER(LEN=15) ::  char_t = 'ls_forcing_top_'   !< leading substring for variables at top boundary 
-    
+
        CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names         !< list of variable in dynamic input file
        CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_l  !< names of mesoscale nested chemistry variables at left boundary
@@ -363,5 +383,5 @@
 
        LOGICAL      ::  init         = .FALSE. !< flag indicating that offline nesting is already initialized
-       
+
        LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_l !< flags inidicating whether left boundary data for chemistry is in dynamic input file  
        LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_n !< flags inidicating whether north boundary data for chemistry is in dynamic input file 
@@ -491,52 +511,52 @@
 
        !-DIMENSIONS
-       INTEGER(iwp)                                 :: nspec=0                   !< number of chem species for which emission values are provided
-       INTEGER(iwp)                                 :: ncat=0                    !< number of emission categories
-       INTEGER(iwp)                                 :: nvoc=0                    !< number of VOCs components
-       INTEGER(iwp)                                 :: npm=0                     !< number of PMs components
-       INTEGER(iwp)                                 :: nnox=2                    !< number of NOx components: NO and NO2
-       INTEGER(iwp)                                 :: nsox=2                    !< number of SOx components: SO and SO4
-       INTEGER(iwp)                                 :: nhoursyear                !< number of hours of a specific year in the HOURLY mode
-                                                                                 !  of the default mode
-       INTEGER(iwp)                                 :: nmonthdayhour             !< number of month days and hours in the MDH mode 
-                                                                                 !  of the default mode
-       INTEGER(iwp)                                 :: dt_emission               !< Number of emissions timesteps for one year 
-                                                                                 !  in the pre-processed emissions case
+       
+       INTEGER(iwp)                                 :: nspec=0            !< no of chem species provided in emission_values
+       INTEGER(iwp)                                 :: n_emiss_species=0  !< no of chem species provided in emission_values
+                                                                          !< same function as nspec, which will be depreciated (ecc)
+                                                                                 
+       INTEGER(iwp)                                 :: ncat=0             !< number of emission categories
+       INTEGER(iwp)                                 :: nvoc=0             !< number of VOC components
+       INTEGER(iwp)                                 :: npm=0              !< number of PM components
+       INTEGER(iwp)                                 :: nnox=2             !< number of NOx components: NO and NO2
+       INTEGER(iwp)                                 :: nsox=2             !< number of SOX components: SO and SO4
+       INTEGER(iwp)                                 :: nhoursyear         !< number of hours of a specific year in the HOURLY mode
+                                                                          !< of the default mode
+       INTEGER(iwp)                                 :: nmonthdayhour      !< number of month days and hours in the MDH mode 
+                                                                          !< of the default mode
+       INTEGER(iwp)                                 :: dt_emission        !< Number of emissions timesteps for one year 
+                                                                          !< in the pre-processed emissions case
        !-- 1d emission input variables
-       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: pm_name                   !< Names of PMs components
-       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: cat_name                  !< Emission categories names
-       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: species_name              !< Names of emission chemical species
-       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: voc_name                  !< Names of VOCs components
-       CHARACTER (LEN=25)                           :: units                     !< Units
-
-       INTEGER(iwp)                                 :: i_hour                    !< indices for assigning the emission values at different timesteps
-       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: cat_index                 !< Index of emission categories
-       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: species_index             !< Index of emission chem species
-
-       REAL(wp),ALLOCATABLE, DIMENSION(:)           :: xm                        !< Molecular masses of emission chem species
+       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: pm_name       !< Names of PM components
+       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: cat_name      !< Emission category names
+       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: species_name  !< Names of emission chemical species
+       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: voc_name      !< Names of VOCs components
+       CHARACTER (LEN=25)                           :: units         !< Units
+
+       INTEGER(iwp)                                 :: i_hour         !< indices for assigning emission values at different timesteps
+       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: cat_index      !< Indices for emission categories
+       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: species_index  !< Indices for emission chem species
+
+       REAL(wp),ALLOCATABLE, DIMENSION(:)           :: xm             !< Molecular masses of emission chem species
 
        !-- 2d emission input variables 
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: hourly_emis_time_factor   !< Time factors for HOURLY emissions (DEFAULT mode)
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: mdh_emis_time_factor      !< Time factors for MDH emissions (DEFAULT mode)
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: nox_comp                  !< Composition of NO and NO2 
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: sox_comp                  !< Composition of SO2 and SO4
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: voc_comp                  !< Composition of different VOC components (number not fixed)
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: hourly_emis_time_factor  !< Time factors for HOURLY emissions (DEFAULT mode)
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: mdh_emis_time_factor     !< Time factors for MDH emissions (DEFAULT mode)
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: nox_comp                 !< Composition of NO and NO2 
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: sox_comp                 !< Composition of SO2 and SO4
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: voc_comp                 !< Composition of VOC components (not fixed)
 
        !-- 3d emission input variables
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)       :: pm_comp                   !< Composition of different PMs components (number not fixed) 
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)       :: pm_comp                  !< Composition of PM components (not fixed) 
   
     END TYPE chem_emis_att_type
 
 
-!-- Data type for the values of chemistry emissions ERUSSO
+!-- Data type for the values of chemistry emissions
     TYPE chem_emis_val_type 
 
-       !REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: stack_height              !< stack height
-
-       !-- 3d emission input variables
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)     :: default_emission_data     !< Input Values emissions DEFAULT mode
-
-       !-- 4d emission input variables 
-       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:)   :: preproc_emission_data      !< Input Values emissions PRE-PROCESSED mode
+       !REAL(wp),ALLOCATABLE, DIMENSION(:,:)     :: stack_height           !< stack height (ecc / to be implemented)
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    :: default_emission_data  !< Emission input values for LOD1 (DEFAULT mode)
+       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:)  :: preproc_emission_data  !< Emission input values for LOD2 (PRE-PROCESSED mode)
 
     END TYPE chem_emis_val_type
@@ -890,5 +910,7 @@
        MODULE PROCEDURE get_variable_4d_real
        MODULE PROCEDURE get_variable_5d_to_4d_real
-       MODULE PROCEDURE get_variable_string       
+       MODULE PROCEDURE get_variable_5d_real           ! (ecc) temp subroutine 4 reading 5D NC arrays
+       MODULE PROCEDURE get_variable_5d_real_dynamic   ! 2B removed as z is out of emission_values
+       MODULE PROCEDURE get_variable_string
     END INTERFACE get_variable
 
@@ -1380,10 +1402,11 @@
 ! Description:
 ! ------------
-!> Reads Chemistry NETCDF Input data, such as emission values, emission species, etc. .
-!------------------------------------------------------------------------------!
+!> Reads Chemistry NETCDF Input data, such as emission values, emission species, etc.
+!------------------------------------------------------------------------------!
+
     SUBROUTINE netcdf_data_input_chemistry_data(emt_att,emt)
 
        USE chem_modules,                                       &
-           ONLY:  mode_emis, time_fac_type, surface_csflux_name
+           ONLY:  emiss_lod, time_fac_type, surface_csflux_name
 
        USE control_parameters,                                 &
@@ -1395,268 +1418,498 @@
        IMPLICIT NONE
 
-       TYPE(chem_emis_att_type), INTENT(INOUT)                                        :: emt_att
-       TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)             :: emt
+       TYPE(chem_emis_att_type), INTENT(INOUT)                             ::  emt_att
+       TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)  ::  emt
     
-       INTEGER(iwp)                                     :: ispec                 !< index for number of emission species in input
-
-       INTEGER(iwp)                                     :: num_vars              !< number of variables in netcdf input file
-       INTEGER(iwp)                                     :: len_dims              !< Length of dimension
-
-       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)          :: dum_var_3d            !< variable for storing temporary data of 3-dimensional 
-                                                                                 !  variables read from netcdf for chemistry emissions
-
-       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)        :: dum_var_4d            !< variable for storing temporary data of 5-dimensional
-                                                                                 !< variables read from netcdf for chemistry emissions
-!--
-       !> Start the processing of the data
-
-       !> Parameterized mode of the emissions
-       IF (TRIM(mode_emis)=="PARAMETERIZED" .OR. TRIM(mode_emis)=="parameterized") THEN 
+       INTEGER(iwp)  ::  i, j, k      !< generic counters
+       INTEGER(iwp)  ::  ispec        !< index for number of emission species in input
+       INTEGER(iwp)  ::  len_dims     !< Length of dimension
+       INTEGER(iwp)  ::  num_vars     !< number of variables in netcdf input file
+
+!
+!-- dum_var_4d are designed to read in emission_values from the chemistry netCDF file.
+!-- Currently the vestigial "z" dimension in emission_values makes it a 5D array,
+!-- hence the corresponding dum_var_5d array.  When the "z" dimension is removed
+!-- completely, dum_var_4d will be used instead
+!-- (ecc 20190425) 
+
+!       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)    ::  dum_var_4d  !< temp array 4 4D chem emission data
+       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:)  ::  dum_var_5d  !< temp array 4 5D chem emission data
+
+!
+!-- Start processing data
+
+       CALL location_message( 'starting allocation of chemistry emissions arrays', .FALSE. )
+
+!
+!-- Emission LOD 0 (Parameterized mode)
+
+        IF  ( emiss_lod == 0 )  THEN
+
+! for reference (ecc)
+!       IF (TRIM(mode_emis) == "PARAMETERIZED" .OR. TRIM(mode_emis) == "parameterized") THEN 
 
            ispec=1
-           emt_att%nspec=0
-
-          !number of species
+           emt_att%n_emiss_species = 0
+
+!
+!-- number of species
+
            DO  WHILE (TRIM( surface_csflux_name( ispec ) ) /= 'novalue' )
 
-             emt_att%nspec=emt_att%nspec+1
+             emt_att%n_emiss_species = emt_att%n_emiss_species + 1
              ispec=ispec+1
+!
+!-- followling line retained for compatibility with salsa_mod
+!-- which still uses emt_att%nspec heavily (ecc)
+
+             emt_att%nspec = emt_att%nspec + 1
 
            ENDDO
 
-          !-- allocate emission values data type arrays
-          ALLOCATE(emt(emt_att%nspec))
-
-          !-- Read EMISSION SPECIES NAMES
-
-          !Assign values
-          ALLOCATE(emt_att%species_name(emt_att%nspec))
+!
+!-- allocate emission values data type arrays
+
+          ALLOCATE ( emt(emt_att%n_emiss_species) )
+
+!
+!-- Read EMISSION SPECIES NAMES
+
+!
+!-- allocate space for strings
+
+          ALLOCATE (emt_att%species_name(emt_att%n_emiss_species) )
  
-         DO ispec=1,emt_att%nspec
+         DO ispec = 1, emt_att%n_emiss_species
             emt_att%species_name(ispec) = TRIM(surface_csflux_name(ispec))
          ENDDO
 
-
-       !> DEFAULT AND PRE-PROCESSED MODE
+!
+!-- LOD 1 (default mode) and LOD 2 (pre-processed mode)
+
        ELSE
 
-#if defined ( __netcdf )        
+#if defined ( __netcdf )
+
           IF ( .NOT. input_pids_chem )  RETURN
 
-          !-- Open file in read-only mode
-          CALL open_read_file( TRIM( input_file_chem ) //                       &
-                               TRIM( coupling_char ), id_emis )
-          !-- inquire number of variables
-          CALL inquire_num_variables( id_emis, num_vars )
-
-          !-- Get General Dimension Lengths: only number of species and number of categories.
-          !                                  the other dimensions depend on the mode of the emissions or on the presence of specific components
-          !nspecies
-          CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nspec, 'nspecies' )
-
+!
+!-- first we allocate memory space for the emission species and then
+!-- we differentiate between LOD 1 (default mode) and LOD 2 (pre-processed mode)
+
+!
+!-- open emission data file ( {palmcase}_chemistry )
+
+          CALL open_read_file ( TRIM(input_file_chem) // TRIM(coupling_char), id_emis )
+
+!
+!-- inquire number of variables
+
+          CALL inquire_num_variables ( id_emis, num_vars )
+
+!
+!-- Get General Dimension Lengths: only # species and # categories.
+!-- Tther dimensions depend on the emission mode or specific components
+
+          CALL netcdf_data_input_get_dimension_length (    &
+                                 id_emis, emt_att%n_emiss_species, 'nspecies' )
+
+!
+!-- backward compatibility for salsa_mod (ecc)
+
+          emt_att%nspec = emt_att%n_emiss_species
+
+! 
+!-- Allocate emission values data type arrays
+
+          ALLOCATE ( emt(emt_att%n_emiss_species) )
+
+!
+!-- READING IN SPECIES NAMES
+
+!
+!-- Allocate memory for species names
+
+          ALLOCATE ( emt_att%species_name(emt_att%n_emiss_species) )
+
+!
+!-- Retrieve variable name (again, should use n_emiss_strlen)
+
+          CALL get_variable( id_emis, 'emission_name',    &
+                             string_values, emt_att%n_emiss_species )
+          emt_att%species_name=string_values
+
+!
+!-- dealocate string_values previously allocated in get_variable call
+
+          IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)
+
+!
+!-- READING IN SPECIES INDICES
+
+!
+!-- Allocate memory for species indices
+
+          ALLOCATE ( emt_att%species_index(emt_att%n_emiss_species) )
+
+!
+!-- Retrieve variable data
+
+          CALL get_variable( id_emis, 'emission_index', emt_att%species_index )
+!
+!-- Now the routine has to distinguish between chemistry emission
+!-- LOD 1 (DEFAULT mode) and LOD 2 (PRE-PROCESSED mode) 
+
+!
+!-- START OF EMISSION LOD 1 (DEFAULT MODE)
+
+
+          IF  ( emiss_lod == 1 )  THEN
+
+! for reference (ecc)
+!          IF (TRIM(mode_emis) == "DEFAULT" .OR. TRIM(mode_emis) == "default") THEN
+
+! 
+!-- get number of emission categories
+
+             CALL netcdf_data_input_get_dimension_length (           &
+                                    id_emis, emt_att%ncat, 'ncat' )
+
+!-- READING IN EMISSION CATEGORIES INDICES
+
+             ALLOCATE ( emt_att%cat_index(emt_att%ncat) )
+
+!
+!-- Retrieve variable data
+
+             CALL get_variable( id_emis, 'emission_cat_index', emt_att%cat_index )
+
+
+!
+!-- Loop through individual species to get basic information on
+!-- VOC/PM/NOX/SOX
+
+!------------------------------------------------------------------------------
+!-- NOTE - CHECK ARRAY INDICES FOR READING IN NAMES AND SPECIES
+!--        IN LOD1 (DEFAULT MODE) FOR THE VARIOUS MODE SPLITS
+!--        AS ALL ID_EMIS CONDITIONALS HAVE BEEN REMOVED FROM GET_VAR
+!--        FUNCTIONS.  IN THEORY THIS WOULD MEAN ALL ARRAYS SHOULD BE
+!--        READ FROM 0 to N-1 (C CONVENTION) AS OPPOSED TO 1 to N
+!--        (FORTRAN CONVENTION).  KEEP THIS IN MIND !!
+!--        (ecc 20190424)
+!------------------------------------------------------------------------------
  
-          !-- Allocate emission values data type arrays
-          ALLOCATE(emt(1:emt_att%nspec))
-
-
-          !-- Read EMISSION SPECIES NAMES
-          !Allocate Arrays 
-          ALLOCATE(emt_att%species_name(emt_att%nspec)) 
-
-          !Call get Variable
-          CALL get_variable( id_emis, 'emission_name', string_values, emt_att%nspec )
-          emt_att%species_name=string_values
-          ! If allocated, Deallocate var_string, an array only used for reading-in strings
-          IF (ALLOCATED(string_values)) DEALLOCATE(string_values) 
-
-          !-- Read EMISSION SPECIES INDEX
-          !Allocate Arrays 
-          ALLOCATE(emt_att%species_index(emt_att%nspec))
-          !Call get Variable
-          CALL get_variable( id_emis, 'emission_index', emt_att%species_index )
-
-
-          !-- Now the routine has to distinguish between DEFAULT and PRE-PROCESSED chemistry emission modes
-
-          IF (TRIM(mode_emis)=="DEFAULT" .OR. TRIM(mode_emis)=="default") THEN
+             DO  ispec = 1, emt_att%n_emiss_species
+
+!
+!-- VOC DATA (name and composition)
+
+                IF  ( TRIM(emt_att%species_name(ispec)) == "VOC" .OR.                  &
+                      TRIM(emt_att%species_name(ispec)) == "voc" )  THEN
+
+!
+!-- VOC name
+                   CALL netcdf_data_input_get_dimension_length (     &
+                                          id_emis, emt_att%nvoc, 'nvoc' )
+                   ALLOCATE ( emt_att%voc_name(emt_att%nvoc) )
+                   CALL get_variable ( id_emis,"emission_voc_name",  &
+                                       string_values, emt_att%nvoc )
+                   emt_att%voc_name = string_values
+                   IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)
+
+! 
+!-- VOC composition
+
+                   ALLOCATE ( emt_att%voc_comp(emt_att%ncat,emt_att%nvoc) )
+                   CALL get_variable ( id_emis, "composition_voc", emt_att%voc_comp,     &
+                                       1, emt_att%ncat, 1, emt_att%nvoc )
+
+                ENDIF  ! VOC
+
+!
+!-- PM DATA (name and composition)
+
+                IF  ( TRIM(emt_att%species_name(ispec)) == "PM" .OR.                   &
+                      TRIM(emt_att%species_name(ispec)) == "pm")  THEN
+
+!
+!-- PM name
+
+                   CALL netcdf_data_input_get_dimension_length (     &
+                                          id_emis, emt_att%npm, 'npm' )
+                   ALLOCATE ( emt_att%pm_name(emt_att%npm) )
+                   CALL get_variable ( id_emis, "pm_name", string_values, emt_att%npm )
+                   emt_att%pm_name = string_values
+                   IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)      
+
+!
+!-- PM composition (PM1, PM2.5 and PM10)
+
+                   len_dims = 3  ! PM1, PM2.5, PM10
+                   ALLOCATE(emt_att%pm_comp(emt_att%ncat,emt_att%npm,len_dims))
+                   CALL get_variable ( id_emis, "composition_pm", emt_att%pm_comp,       &
+                                       1, emt_att%ncat, 1, emt_att%npm, 1, len_dims )
+
+                ENDIF  ! PM
+
+!
+!-- NOX (NO and NO2)
+
+                IF  ( TRIM(emt_att%species_name(ispec)) == "NOX" .OR.                  &
+                      TRIM(emt_att%species_name(ispec)) == "nox" )  THEN
+
+                   ALLOCATE ( emt_att%nox_comp(emt_att%ncat,emt_att%nnox) )
+                   CALL get_variable ( id_emis, "composition_nox", emt_att%nox_comp,     &
+                                       1, emt_att%ncat, 1, emt_att%nnox )
+
+                ENDIF  ! NOX
+
+!
+!-- SOX (SO2 and SO4)
+
+                IF  ( TRIM(emt_att%species_name(ispec)) == "SOX" .OR.                  &
+                      TRIM(emt_att%species_name(ispec)) == "sox" )  THEN
+
+                   ALLOCATE ( emt_att%sox_comp(emt_att%ncat,emt_att%nsox) )
+                   CALL get_variable ( id_emis, "composition_sox", emt_att%sox_comp,     &
+                                       1, emt_att%ncat, 1, emt_att%nsox )
+
+                ENDIF  ! SOX
+
+             ENDDO  ! do ispec
+
+!
+!-- EMISSION TIME SCALING FACTORS (hourly and MDH data)
  
-             !number of categories
-             CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%ncat, 'ncat' )
-
-             !-- Read EMISSION CATEGORIES INDEX
-             !Allocate Arrays 
-             ALLOCATE(emt_att%cat_index(emt_att%ncat))
-             !Call get Variable
-             CALL get_variable( id_emis, 'emission_cat_index', emt_att%cat_index )
-
+!     
+!-- HOUR   
+             IF  ( TRIM(time_fac_type) == "HOUR" .OR.                        &
+                   TRIM(time_fac_type) == "hour" )  THEN
+
+                CALL netcdf_data_input_get_dimension_length (                  &
+                                       id_emis, emt_att%nhoursyear, 'nhoursyear' )
+                ALLOCATE ( emt_att%hourly_emis_time_factor(emt_att%ncat,emt_att%nhoursyear) )
+                CALL get_variable ( id_emis, "emission_time_factors",          &
+                                    emt_att%hourly_emis_time_factor,           &
+                                    1, emt_att%ncat, 1, emt_att%nhoursyear )
+
+!
+!-- MDH
+
+             ELSE IF  ( TRIM(time_fac_type)  ==  "MDH" .OR.                  &
+                        TRIM(time_fac_type)  ==  "mdh" )  THEN
+
+                CALL netcdf_data_input_get_dimension_length (                  &
+                                       id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' )
+                ALLOCATE ( emt_att%mdh_emis_time_factor(emt_att%ncat,emt_att%nmonthdayhour) )
+                CALL get_variable ( id_emis, "emission_time_factors",          &
+                                    emt_att%mdh_emis_time_factor,              &
+                                    1, emt_att%ncat, 1, emt_att%nmonthdayhour )
+
+!
+!-- ERROR (time factor undefined)
+
+             ELSE
+
+                message_string = 'We are in the DEFAULT chemistry emissions mode: '  //  &
+                                 '     !no time-factor type specified!'              //  &
+                                 'Please specify the value of time_fac_type:'        //  &
+                                 '         either "MDH" or "HOUR"'                  
+                CALL message( 'netcdf_data_input_chemistry_data', 'CM0200', 2, 2, 0, 6, 0 ) 
  
-             DO ispec=1,emt_att%nspec
-                !-- EMISSION_VOC_NAME (1-DIMENSIONAL)
-                IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN
-                   !Allocate Array
-                   CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
-                   ALLOCATE(emt_att%voc_name(1:emt_att%nvoc))
-                   !Read-in Variable
-                   CALL get_variable( id_emis,"emission_voc_name",string_values, emt_att%nvoc)
-                   emt_att%voc_name=string_values
-                   IF (ALLOCATED(string_values)) DEALLOCATE(string_values)
+
+             ENDIF  ! time_fac_type
+
+!
+!-- read in default (LOD1) emissions from chemisty netCDF file per species
+
+!
+!-- NOTE - at the moment the data is read in per species, but in the future it would
+!--        be much more sensible to read in per species per time step to reduce
+!--        memory consumption and, to a lesser degree, dimensionality of data exchange
+!--        (I expect this will be necessary when the problem size is large)
+
+             DO ispec = 1, emt_att%n_emiss_species
+
+!
+!-- allocate space for species specific emission values
+!-- NOTE - this array is extended by 1 cell in each horizontal direction
+!--        to compensate for an apparent linear offset.  The reason of this
+!--        offset is not known but it has been determined to take place beyond the
+!--        scope of this module, and has little to do with index conventions.
+!--        That is, setting the array horizontal limit from nx0:nx1 to 1:(nx1-nx0+1)
+!--        or nx0+1:nx1+1 did not result in correct or definite behavior
+!--        This must be looked at at some point by the Hannover team but for now
+!--        this workaround is deemed reasonable (ecc 20190417)
+
+                IF ( .NOT. ALLOCATED ( emt(ispec)%default_emission_data ) )  THEN
+                    ALLOCATE ( emt(ispec)%default_emission_data(emt_att%ncat,nys:nyn+1,nxl:nxr+1) )
+                ENDIF
+!
+!-- allocate dummy variable w/ index order identical to that shown in the netCDF header
+
+                ALLOCATE ( dum_var_5d(1,nys:nyn,nxl:nxr,1,emt_att%ncat) )
+!
+!-- get variable.  be very careful
+!-- I am using get_variable_5d_real_dynamic (note logical argument at the end)
+!-- 1) use Fortran index convention (i.e., 1 to N)
+!-- 2) index order must be in reverse order from above allocation order
  
-                !-- COMPOSITION VOC (2-DIMENSIONAL)
-                   !Allocate Array 
-                   ALLOCATE(emt_att%voc_comp(1:emt_att%ncat,1:emt_att%nvoc))
-                   !Read-in Variable
-!               CALL get_variable(id_emis,"composition_voc",emt%voc_comp,1,1,emt%ncat,emt%nvoc)
-                   CALL get_variable(id_emis,"composition_voc",emt_att%voc_comp,1,emt_att%ncat,1,emt_att%nvoc)
+                CALL get_variable ( id_emis, "emission_values", dum_var_5d, &
+                                    1,            ispec, nxl+1,     nys+1,     1,                    &
+                                    emt_att%ncat, 1,     nxr-nxl+1, nyn-nys+1, emt_att%dt_emission,  &
+                                    .FALSE. )
+!
+!-- assign temp array to data structure then deallocate temp array
+!-- NOTE - indices are shifted from nx0:nx1 to nx0+1:nx1+1 to offset
+!--        the emission data array to counter said domain offset
+!--        (ecc 20190417)
+
+                DO k = 1, emt_att%ncat
+                   DO j = nys+1, nyn+1
+                      DO i = nxl+1, nxr+1
+                         emt(ispec)%default_emission_data(k,j,i) = dum_var_5d(1,j-1,i-1,1,k)
+                      ENDDO
+                   ENDDO
+                ENDDO
+
+                DEALLOCATE ( dum_var_5d )
+
+             ENDDO  ! ispec
+!
+!-- UNITS
+
+             CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values")
+
+!
+!-- END DEFAULT MODE
+
+
+!
+!-- START LOD 2 (PRE-PROCESSED MODE)
+
+          ELSE IF  ( emiss_lod == 2 )  THEN
+
+! for reference (ecc)
+!          ELSE IF (TRIM(mode_emis) == "PRE-PROCESSED" .OR. TRIM(mode_emis) == "pre-processed") THEN
+
+!
+!-- For LOD 2 only VOC and emission data need be read
+
+!------------------------------------------------------------------------------
+!-- NOTE - CHECK ARRAY INDICES FOR READING IN NAMES AND SPECIES
+!--        IN LOD2 (PRE-PROCESSED MODE) FOR THE VARIOUS MODE SPLITS
+!--        AS ALL ID_EMIS CONDITIONALS HAVE BEEN REMOVED FROM GET_VAR
+!--        FUNCTIONS.  IN THEORY THIS WOULD MEAN ALL ARRAYS SHOULD BE
+!--        READ FROM 0 to N-1 (C CONVENTION) AS OPPOSED TO 1 to N
+!--        (FORTRAN CONVENTION).  KEEP THIS IN MIND !!
+!--        (ecc 20190424)
+!------------------------------------------------------------------------------
+
+             DO ispec = 1, emt_att%n_emiss_species
+
+!
+!-- VOC DATA (name and composition)
+
+                IF  ( TRIM(emt_att%species_name(ispec)) == "VOC" .OR.                  &
+                      TRIM(emt_att%species_name(ispec)) == "voc" )  THEN
+
+!
+!-- VOC name
+                   CALL netcdf_data_input_get_dimension_length (                         &
+                                          id_emis, emt_att%nvoc, 'nvoc' )
+                   ALLOCATE ( emt_att%voc_name(emt_att%nvoc) )
+                   CALL get_variable ( id_emis, "emission_voc_name",                     &
+                                       string_values, emt_att%nvoc)
+                   emt_att%voc_name = string_values
+                   IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)
+
+!
+!-- VOC composition
+ 
+                   ALLOCATE ( emt_att%voc_comp(emt_att%ncat,emt_att%nvoc) )
+                   CALL get_variable ( id_emis, "composition_voc", emt_att%voc_comp,     &
+                                       1, emt_att%ncat, 1, emt_att%nvoc )
+                ENDIF  ! VOC
+  
+             ENDDO  ! ispec
+
+!
+!-- EMISSION DATA
+
+             CALL netcdf_data_input_get_dimension_length (                               &
+                                    id_emis, emt_att%dt_emission, 'time' )   
+ 
+!
+!-- read in pre-processed (LOD2) emissions from chemisty netCDF file per species
+
+!
+!-- NOTE - at the moment the data is read in per species, but in the future it would
+!--        be much more sensible to read in per species per time step to reduce
+!--        memory consumption and, to a lesser degree, dimensionality of data exchange
+!--        (I expect this will be necessary when the problem size is large)
+
+             DO ispec = 1, emt_att%n_emiss_species
+
+!
+!-- allocate space for species specific emission values
+!-- NOTE - this array is extended by 1 cell in each horizontal direction
+!--        to compensate for an apparent linear offset.  The reason of this
+!--        offset is not known but it has been determined to take place beyond the
+!--        scope of this module, and has little to do with index conventions.
+!--        That is, setting the array horizontal limit from nx0:nx1 to 1:(nx1-nx0+1)
+!--        or nx0+1:nx1+1 did not result in correct or definite behavior
+!--        This must be looked at at some point by the Hannover team but for now
+!--        this workaround is deemed reasonable (ecc 20190417)
+
+                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) )  THEN
+                   ALLOCATE( emt(ispec)%preproc_emission_data(                           &
+                             emt_att%dt_emission, 1, nys:nyn+1, nxl:nxr+1) )
                 ENDIF
-
-                !-- EMISSION_PM_NAME (1-DIMENSIONAL)
-                IF (TRIM(emt_att%species_name(ispec))=="PM" .OR. TRIM(emt_att%species_name(ispec))=="pm") THEN
-                   CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%npm, 'npm' )
-                   ALLOCATE(emt_att%pm_name(1:emt_att%npm))
-                   !Read-in Variable
-                   CALL get_variable( id_emis,"pm_name",string_values, emt_att%npm)
-                   emt_att%pm_name=string_values
-                   IF (ALLOCATED(string_values)) DEALLOCATE(string_values)      
-
-                !-- COMPOSITION PM (3-DIMENSIONAL)
-                   !Allocate
-                   len_dims=3  !> number of PMs: PM1, PM2.5 and PM10
-                   ALLOCATE(emt_att%pm_comp(1:emt_att%ncat,1:emt_att%npm,1:len_dims))
-                   !Read-in Variable
-                   CALL get_variable(id_emis,"composition_pm",emt_att%pm_comp,1,emt_att%ncat,1,emt_att%npm,1,len_dims)                    
-                ENDIF
-
-                !-- COMPOSITION_NOX (2-DIMENSIONAL)
-                IF (TRIM(emt_att%species_name(ispec))=="NOx" .OR. TRIM(emt_att%species_name(ispec))=="nox") THEN
-                   !Allocate array
-                   ALLOCATE(emt_att%nox_comp(1:emt_att%ncat,1:emt_att%nnox))
-                   !Read-in Variable
-                   CALL get_variable(id_emis,"composition_nox",emt_att%nox_comp,1,emt_att%ncat,1,emt_att%nnox)
-                ENDIF
-
-                !-- COMPOSITION-SOX (2-DIMENSIONAL)
-                IF (TRIM(emt_att%species_name(ispec))=="SOx" .OR. TRIM(emt_att%species_name(ispec))=="sox") THEN
-                   ALLOCATE(emt_att%sox_comp(1:emt_att%ncat,1:emt_att%nsox))
-                   !Read-in Variable
-                   CALL get_variable(id_emis,"composition_sox",emt_att%sox_comp,1,emt_att%ncat,1,emt_att%nsox)
-                ENDIF
-             ENDDO !>ispec
-
-!-- For reading the emission time factors, the distinction between HOUR and MDH data is necessary
-      
-             !-- EMISSION_TIME_SCALING_FACTORS
-                !-- HOUR   
-             IF (TRIM(time_fac_type)=="HOUR" .OR. TRIM(time_fac_type)=="hour") THEN
-                !-- Allocate Array
-                CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nhoursyear, 'nhoursyear' )
-                ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear))
-                !Read-in Variable
-                CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,   &
-                                  emt_att%ncat,1,emt_att%nhoursyear)
-
-                !-- MDH
-             ELSE IF (TRIM(time_fac_type) == "MDH" .OR. TRIM(time_fac_type) == "mdh") THEN
-                !-- Allocate Array
-                CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' )
-                ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour))
-                !-- Read-in Variable
-                CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,       &
-                                  emt_att%ncat,1,emt_att%nmonthdayhour)
-
-             ELSE
-
-             message_string = 'We are in the DEFAULT chemistry emissions mode: '            //          &
-                              '     !no time-factor type specified!'                        //          &
-                              'Please, specify the value of time_fac_type:'                 //          &
-                              '         either "MDH" or "HOUR"'                  
-             CALL message( 'netcdf_data_input_chemistry_data', 'CM0200', 2, 2, 0, 6, 0 ) 
- 
-
-             ENDIF
-
-             !-- Finally read-in the emission values and their units (DEFAULT mode)
-
-             DO ispec=1,emt_att%nspec
-
-                IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) )                              &
-                    ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1))
-
-                ALLOCATE(dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1))
-
-                CALL get_variable(id_emis,"emission_values",dum_var_3d,ispec,1,emt_att%ncat,nys,nyn,nxl,nxr)          
-
-                emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1) =                           &
-                    dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)
-
-                DEALLOCATE (dum_var_3d)
-
-             ENDDO
-
-             !-- UNITS
-             CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values")
-
-
-          !-- PRE-PROCESSED MODE --
-
-          ELSE IF (TRIM(mode_emis)=="PRE-PROCESSED" .OR. TRIM(mode_emis)=="pre-processed") THEN
-          !-- In the PRE-PROCESSED mode, only the VOC names, the VOC_composition, the emission values and their units remain to be read at this point
-
-             DO ispec=1,emt_att%nspec
-
-             !-- EMISSION_VOC_NAME (1-DIMENSIONAL)
-                IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN
-                   !Allocate Array
-                   CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
-                   ALLOCATE(emt_att%voc_name(1:emt_att%nvoc))
-                   !Read-in Variable
-                   CALL get_variable( id_emis,"emission_voc_name",string_values, emt_att%nvoc)
-                   emt_att%voc_name=string_values
-                   IF (ALLOCATED(string_values)) DEALLOCATE(string_values)
- 
-             !-- COMPOSITION VOC (2-DIMENSIONAL)
-                   !Allocate Array 
-                   ALLOCATE(emt_att%voc_comp(1:emt_att%ncat,1:emt_att%nvoc))
-                   !Read-in Variable
-                   CALL get_variable(id_emis,"composition_voc",emt_att%voc_comp,1,emt_att%ncat,1,emt_att%nvoc)
-                ENDIF
-  
-             ENDDO !> ispec
-
-             !-- EMISSION_VALUES (4-DIMENSIONAL)
-             !Calculate temporal dimension length
-             CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%dt_emission, 'time' )   
-         
-
-             DO ispec=1,emt_att%nspec
-
-                !Allocation for the entire domain has to be done only for the first processor between all the subdomains     
-                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) )                              &
-                    ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1))
-
-                !> allocate variable where to pass emission values read from netcdf
-                ALLOCATE(dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1))
-
-                !Read-in Variable 
-                CALL get_variable(id_emis,"emission_values",dum_var_4d,ispec,1,emt_att%dt_emission,1,1,nys,nyn,nxl,nxr)          
-
-      
-                emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1) =                         &
-                      dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1)
-
-                DEALLOCATE ( dum_var_4d )
-
-             ENDDO
-
-             !-- UNITS
-             CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values")
+!
+!-- allocate dummy variable w/ index order identical to that shown in the netCDF header
+
+                ALLOCATE ( dum_var_5d(emt_att%dt_emission,1,nys:nyn,nxl:nxr,1) )
+!
+!-- get variable.  be very careful
+!-- I am using get_variable_5d_real_dynamic (note logical argument at the end)
+!-- 1) use Fortran index convention (i.e., 1 to N)
+!-- 2) index order must be in reverse order from above allocation order
+
+                CALL get_variable ( id_emis, "emission_values", dum_var_5d, &
+                                    ispec, nxl+1,     nys+1,     1, 1,                   &
+                                    1,     nxr-nxl+1, nyn-nys+1, 1, emt_att%dt_emission, &
+                                    .FALSE. )
+!
+!-- assign temp array to data structure then deallocate temp array
+!-- NOTE - indices are shifted from nx0:nx1 to nx0+1:nx1+1 to offset
+!--        the emission data array to counter said unkonwn offset
+!--        (ecc 20190417)
+
+                DO k = 1, emt_att%dt_emission
+                   DO j = nys+1, nyn+1
+                      DO i = nxl+1, nxr+1
+                         emt(ispec)%preproc_emission_data(k,1,j,i) = dum_var_5d(k,1,j-1,i-1,1)
+                      ENDDO
+                   ENDDO
+                ENDDO
+
+                DEALLOCATE ( dum_var_5d )
+
+             ENDDO  ! ispec
+!
+!-- UNITS
+
+             CALL get_attribute ( id_emis, "units", emt_att%units, .FALSE. , "emission_values" )
        
-          ENDIF
-
-       CALL close_input_file( id_emis )
+          ENDIF  ! LOD1 & LOD2 (default and pre-processed mode)
+
+          CALL close_input_file (id_emis)
 
 #endif
-       ENDIF
+
+       ENDIF ! LOD0 (parameterized mode)
 
     END SUBROUTINE netcdf_data_input_chemistry_data
+
 
 !------------------------------------------------------------------------------!
@@ -5358,36 +5611,15 @@
        ENDIF
 
-
-       !Temporary solution for reading emission chemistry files: TBD: we should discuss whether remove it or not
-       IF ( id==id_emis ) THEN
-
-          !--    Allocate temporary variable according to memory access on file. 
-          ALLOCATE( tmp(is:ie,js:je) )
-
-          !--    Get variable
-          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
-                                  start = (/ is,      js /),                  &
-                                  count = (/ ie-is+1 , je-js+1 /) ) 
-
-          var=tmp
-
-          CALL handle_error( 'get_variable_2d_real', 530, variable_name ) !TBD: the error number shuld be changed, but since the solution is 
-                                                                          ! provisory, we give the same as below
- 
-          DEALLOCATE( tmp )
-       
-       !>  Original Subroutine part
-       ELSE
-!
-!--       Allocate temporary variable according to memory access on file. 
-          ALLOCATE( tmp(is:ie,js:je) )
-!
-!--       Get variable
-          nc_stat = NF90_GET_VAR( id, id_var, tmp,                             &
-                                  start = (/ is+1,      js+1 /),               &
-                                  count = (/ ie-is + 1, je-js+1 /) )   
+!
+!-- Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je) )
+!
+!-- Get variable
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,            &
+                      start = (/ is+1,      js+1 /),       &
+                      count = (/ ie-is + 1, je-js+1 /) )   
           CALL handle_error( 'get_variable_2d_real', 530, variable_name )
 !
-!--       Resort data. Please note, dimension subscripts of var all start at 1. 
+!-- Resort data. Please note, dimension subscripts of var all start at 1. 
           DO  i = is, ie 
              DO  j = js, je 
@@ -5398,5 +5630,4 @@
           DEALLOCATE( tmp )
 
-       ENDIF
 #endif
     END SUBROUTINE get_variable_2d_real
@@ -5713,49 +5944,30 @@
        ENDIF
 
-      !Temporary solution for reading emission chemistry files: 
-       IF ( id==id_emis ) THEN
-
-          !--    Allocate temporary variable according to memory access on file. 
-          ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) )
-
-          !--    Get variable
-          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
-                                  start = (/ is,   js,   k1s+1,   k2s+1 /),                  &
-                                  count = (/ ie-is+1 , je-js+1, k1e-k1s+1, k2e-k2s+1 /) ) 
-
-          var=tmp
+!
+!-- Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) )
+!
+!-- Get variable
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                      start = (/ is+1,    js+1,    k1s+1,     k2s+1 /),        &
+                      count = (/ ie-is+1, je-js+1, k1e-k1s+1, k2e-k2s+1 /) )
 
           CALL handle_error( 'get_variable_4d_real', 535, variable_name )
- 
-          DEALLOCATE( tmp )
-
-       !> Original subroutine part
-       ELSE
-!
-!--       Allocate temporary variable according to memory access on file. 
-          ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) )
-!
-!--       Get variable
-          nc_stat = NF90_GET_VAR( id, id_var, tmp,                             &
-                               start = (/ is+1,    js+1,   k1s+1, k2s+1 /),    &
-                               count = (/ ie-is+1, je-js+1,                    &
-                                          k1e-k1s+1, k2e-k2s+1 /) )
-
-          CALL handle_error( 'get_variable_4d_real', 535, variable_name )
-!
-!--       Resort data. Please note, dimension subscripts of var all start at 1. 
-          DO  i = is, ie 
-             DO  j = js, je
-                DO  k1 = k1s, k1e
-                   DO  k2 = k2s, k2e
-                      var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2)
-                   ENDDO
+!
+!-- Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je
+             DO  k1 = k1s, k1e
+                DO  k2 = k2s, k2e
+                   var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2)
                 ENDDO
              ENDDO
           ENDDO
-
-          DEALLOCATE( tmp )
-       ENDIF
+       ENDDO
+
+       DEALLOCATE( tmp )
+
 #endif
+
     END SUBROUTINE get_variable_4d_real
 
@@ -5806,5 +6018,5 @@
 
       !Temporary solution for reading emission chemistry files: 
-       IF ( id==id_emis ) THEN
+       IF ( id == id_emis ) THEN
 
           !--    Allocate temporary variable according to memory access on file. 
@@ -5993,5 +6205,5 @@
 
       !Temporary solution for reading emission chemistry files: 
-       IF ( id==id_emis ) THEN
+       IF ( id == id_emis ) THEN
 
           !--    Allocate temporary variable according to memory access on file. 
@@ -6040,4 +6252,232 @@
     END SUBROUTINE get_variable_5d_to_4d_real
 
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads a 5D float variable from file. 
+!> NOTE - This subroutine is used specific for reading NC variable
+!>        emission_values having a "z" dimension.  Said dimension
+!>        is to be removed in the future and this subroutine shall
+!>        be depreciated accordingly (ecc 20190418)
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_variable_5d_real( id, variable_name, var, is, ie, js, je,   &
+                                     k1s, k1e, k2s, k2e, k3s, k3e )
+
+       USE indices
+       USE pegrid
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*)          ::  variable_name  !< variable name
+
+       INTEGER(iwp)              :: i       !< i index
+       INTEGER(iwp)              :: ie      !< i index start
+       INTEGER(iwp)              :: is      !< i index end
+       INTEGER(iwp)              :: id_var  !< netCDF variable ID (varid)
+       INTEGER(iwp)              :: j       !< j index
+       INTEGER(iwp)              :: je      !< j index start
+       INTEGER(iwp)              :: js      !< j index end
+       INTEGER(iwp)              :: k1      !< k1 index
+       INTEGER(iwp)              :: k1e     !< k1 index start
+       INTEGER(iwp)              :: k1s     !< k1 index end
+       INTEGER(iwp)              :: k2      !< k2 index
+       INTEGER(iwp)              :: k2e     !< k2 index start
+       INTEGER(iwp)              :: k2s     !< k2 index end
+       INTEGER(iwp)              :: k3      !< k3 index
+       INTEGER(iwp)              :: k3e     !< k3 index start
+       INTEGER(iwp)              :: k3s     !< k3 index end
+       INTEGER(iwp), INTENT(IN)  :: id      !< netCDF file ID (ncid)
+
+       REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE    :: tmp  !< temp array to read data from file
+       REAL(wp), DIMENSION(:,:,:,:,:), INTENT(INOUT)  :: var  !< variable to be read 
+
+#if defined( __netcdf )
+
+!
+!-- Inquire variable id
+
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
+
+!
+!-- Check for collective read-operation and set respective NetCDF flags if required.
+ 
+       IF ( collective_read )  THEN
+
+#if defined( __netcdf4_parallel )       
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+#endif
+
+       ENDIF
+
+!
+!-- Allocate temporary variable according to memory access on file. 
+
+       ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e,k3s:k3e) )
+
+!
+!-- Get variable from file
+
+       nc_stat = NF90_GET_VAR ( id, id_var, tmp,                                         &
+                      start = (/ is+1,    js+1,    k1s+1,     k2s+1,     k3s+1 /),       &
+                      count = (/ ie-is+1, je-js+1, k1e-k1s+1, k2e-k2s+1, k3e-k3s+1 /) )
+
+       CALL handle_error( 'get_variable_5d_real', 535, variable_name )
+
+!
+!-- Resort (reverse index order) and standardize (from 1 to N) output array
+
+       DO  i = is, ie 
+          DO  j = js, je
+             DO  k1 = k1s, k1e
+                DO  k2 = k2s, k2e
+                   DO k3 = k3s, k3e
+                      var(k3-k3s+1,k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2,k3)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDDO
+
+       DEALLOCATE( tmp )
+
+#endif
+
+    END SUBROUTINE get_variable_5d_real
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads a 5D float variables from dynamic driver, such as time-dependent xy-, 
+!> xz- or yz-boundary data as well as 5D initialization data. Please note, 
+!> the passed arguments are start indices and number of elements in each 
+!> dimension, which is in contrast to the other 3d versions where start- and 
+!> end indices are passed. The different handling of 5D dynamic variables is 
+!> due to its asymmetry for the u- and v component.
+!> NOTE(1) - This subroutine is more flexible than get_variable_xd_real as it
+!>           provides much better control over starting and count indices
+!>           (ecc 20190418)
+!> NOTE(2) - This subroutine is used specific for reading NC variable
+!>           emission_values having a "z" dimension.  Said dimension
+!>           is to be removed in the future and this subroutine shall
+!>           be depreciated accordingly (ecc 20190418)
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE get_variable_5d_real_dynamic( id, variable_name, var,                       &
+                                             i1s, i2s, i3s, i4s, i5s,                      &
+                                             count_1, count_2, count_3, count_4, count_5,  &
+                                             par_access )
+
+       USE indices
+       USE pegrid
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*)          ::  variable_name  !< variable name
+
+       LOGICAL                   ::  par_access     !< additional flag indicating parallel read
+
+       INTEGER(iwp)              ::  count_1  !< # elements read in dimension 1 wrt file
+       INTEGER(iwp)              ::  count_2  !< # elements read in dimension 2 wrt file
+       INTEGER(iwp)              ::  count_3  !< # elements read in dimension 3 wrt file
+       INTEGER(iwp)              ::  count_4  !< # elements read in dimension 4 wrt file
+       INTEGER(iwp)              ::  count_5  !< # elements read in dimension 5 wrt file
+       INTEGER(iwp)              ::  i1       !< index for dimension 1 on file
+       INTEGER(iwp)              ::  i1s      !< starting index for dimension 1 hyperslab
+       INTEGER(iwp)              ::  i2       !< index for dimension 2 on file
+       INTEGER(iwp)              ::  i2s      !< starting index for dimension 2 hyperslab
+       INTEGER(iwp)              ::  i3       !< index for dimension 3 on file
+       INTEGER(iwp)              ::  i3s      !< starting index for dimension 3 hyperslab
+       INTEGER(iwp)              ::  i4       !< index for dimension 4 on file
+       INTEGER(iwp)              ::  i4s      !< starting index for dimension 4 hyperslab
+       INTEGER(iwp)              ::  i5       !< index for dimension 5 on file
+       INTEGER(iwp)              ::  i5s      !< starting index for dimension 5 hyperslab
+       INTEGER(iwp)              ::  id_var   !< netCDF variable id (varid)
+       INTEGER(iwp)              ::  lb1      !< lower bound of dimension 1 wrt file
+       INTEGER(iwp)              ::  lb2      !< lower bound of dimension 2 wrt file
+       INTEGER(iwp)              ::  lb3      !< lower bound of dimension 3 wrt file
+       INTEGER(iwp)              ::  lb4      !< lower bound of dimension 4 wrt file
+       INTEGER(iwp)              ::  lb5      !< lower bound of dimension 5 wrt file
+       INTEGER(iwp)              ::  ub1      !< upper bound of dimension 1 wrt file
+       INTEGER(iwp)              ::  ub2      !< upper bound of dimension 2 wrt file
+       INTEGER(iwp)              ::  ub3      !< upper bound of dimension 3 wrt file
+       INTEGER(iwp)              ::  ub4      !< upper bound of dimension 4 wrt file
+       INTEGER(iwp)              ::  ub5      !< upper bound of dimension 5 wrt file
+       INTEGER(iwp), INTENT(IN)  ::  id       !< netCDF file id (ncid)
+
+       REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE    ::  tmp  !< temporary variable to read data
+                                                               !< from file according is reverse
+                                                               !< array index order
+       REAL(wp), DIMENSION(:,:,:,:,:), INTENT(INOUT)  ::  var  !< input variable
+       
+#if defined( __netcdf )
+
+!
+!-- Inquire variable id.
+
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
+
+!
+!-- Check for collective read-operation and set respective NetCDF flags if required. 
+!-- Please note, in contrast to the other input routines where each PEs 
+!-- reads its subdomain data, dynamic input data not by all PEs, only
+!-- by those which encompass lateral model boundaries. Hence, collective
+!-- read operations are only enabled for top-boundary data.
+
+       IF ( collective_read  .AND.  par_access )  THEN
+
+#if defined( __netcdf4_parallel )       
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+#endif
+
+       ENDIF
+
+!
+!-- Allocate temporary variable according to memory access on file. 
+!-- Therefore, determine dimension bounds of input array. 
+
+       lb1 = LBOUND(var,5)
+       ub1 = UBOUND(var,5)
+       lb2 = LBOUND(var,4)
+       ub2 = UBOUND(var,4)
+       lb3 = LBOUND(var,3)
+       ub3 = UBOUND(var,3)
+       lb4 = LBOUND(var,2)
+       ub4 = UBOUND(var,2)
+       lb5 = LBOUND(var,1)
+       ub5 = UBOUND(var,1)
+       ALLOCATE ( tmp(lb1:ub1,lb2:ub2,lb3:ub3,lb4:ub4,lb5:ub5) )
+
+!
+!-- Get variable
+
+       nc_stat = NF90_GET_VAR(  id, id_var, tmp,                                         &
+                      start = (/ i1s,     i2s,     i3s,     i4s,     i5s     /),         &
+                      count = (/ count_1, count_2, count_3, count_4, count_5 /) )
+
+       CALL handle_error( 'get_variable_3d_real_dynamic', 537, variable_name )
+
+!
+!-- Assign temp array to output.  Note reverse index order 
+
+       DO  i5 = lb5, ub5
+          DO  i4 = lb4, ub4
+             DO  i3 = lb3, ub3
+                DO i2 = lb2, ub2
+                   DO  i1 = lb1, ub1
+                      var(i5,i4,i3,i2,i1) = tmp(i1,i2,i3,i4,i5)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDDO
+
+       DEALLOCATE( tmp )
+
+#endif
+
+    END SUBROUTINE get_variable_5d_real_dynamic
+
 
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3967)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3968)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! replace nspec_out with n_matched_vars
+! 
+! 3929 2019-04-24 12:52:08Z banzhafs
 ! Reverse changes back from revision 3878: use chem_boundary_conds instead of
 ! chem_boundary_conds_decycle
@@ -518,5 +521,5 @@
 
     USE chem_modules,                                                                              &
-        ONLY:  bc_cs_t_val, chem_species, cs_name, emissions_anthropogenic, nspec_out
+        ONLY:  bc_cs_t_val, chem_species, cs_name, emissions_anthropogenic, n_matched_vars
 
     USE chemistry_model_mod,                                                                       &
@@ -1448,5 +1451,5 @@
 !--       Call emission routine only once an hour
           IF (hour_of_year  .GT.  hour_call_emis )  THEN
-             CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
+             CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
              hour_call_emis = hour_of_year
           ENDIF
