Index: palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3610)
+++ palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3611)
@@ -22,13 +22,19 @@
 ! Current revisions:
 ! ------------------
-! 
+!
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Code update to comply PALM coding rules
+! removed unnecessary informative messsages/location
+! messages and corrected comments on PM units from to kg  
+! bug fix: spcs_name replaced by nvar in DO loops
+!
+! 3591 2018-11-30 17:37:32Z suehring
 ! - Removed salsa dependency.
 ! - Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but 
 !   salsa is (M. Kurppa)
-! 
+!
 ! 3582 2018-11-29 19:16:36Z suehring
 ! resler:
@@ -67,5 +73,5 @@
 ! Authors:
 ! --------
-! @author Emmanuele Russo (Fu-Berlin)
+! @author Emmanuele Russo (FU-Berlin)
 ! @author Sabine Banzhaf  (FU-Berlin)
 ! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
@@ -79,4 +85,6 @@
 !> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
 !> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
+!> @todo fill/activate restart option in chem_emissions_init
+!> @todo discuss dt_emis
 !> @note <Enter notes on the module>
 !> @bug  <Enter known bugs here>
@@ -86,9 +94,9 @@
 
     USE arrays_3d,                                                             &
-       ONLY: rho_air
+        ONLY:  rho_air
 
     USE control_parameters,                                                    &
-       ONLY:  initializing_actions, end_time, message_string,                  &
-              intermediate_timestep_count, dt_3d
+        ONLY:  end_time, message_string, initializing_actions,                 &
+               intermediate_timestep_count, dt_3d
  
     USE indices
@@ -96,148 +104,72 @@
     USE kinds
 
-#if defined ( __netcdf )
-    USE NETCDF
+#if defined( __netcdf )
+    USE netcdf
 #endif
 
     USE netcdf_data_input_mod,                                                  &
-       ONLY: chem_emis_att_type, chem_emis_val_type
+        ONLY: chem_emis_att_type, chem_emis_val_type
 
     USE date_and_time_mod,                                                      &
-       ONLY: time_default_indices, time_preprocessed_indices,                  &
-             month_of_year, day_of_month, hour_of_day,                          &
-             index_mm, index_dd, index_hh
-
+        ONLY: day_of_month, hour_of_day,                                        &
+             index_mm, index_dd, index_hh,                                      &
+             month_of_year, hour_of_day,                                        &
+             time_default_indices, time_preprocessed_indices
+    
     USE chem_gasphase_mod,                                                      &
-       ONLY: spc_names
+        ONLY: spc_names
  
     USE chem_modules
 
     USE statistics,                                                             &   
-           ONLY:  weight_pres
-
+        ONLY:  weight_pres
+
+    
     IMPLICIT NONE
 
+!
 !-- Declare all global variables within the module 
     
-    CHARACTER (LEN=80)                               :: filename_emis                   !> Variable for the name of the netcdf input file
-
-    INTEGER(iwp)                                     :: i                               !> index 1st selected dimension (some dims are not spatial)
-    INTEGER(iwp)                                     :: j                               !> index 2nd selected dimension 
-    INTEGER(iwp)                                     :: i_start                         !> Index to start read variable from netcdf along one dims
-    INTEGER(iwp)                                     :: i_end                           !> Index to end read variable from netcdf in one dims
-    INTEGER(iwp)                                     :: j_start                         !> Index to start read variable from netcdf in additional dims
-    INTEGER(iwp)                                     :: j_end                           !> Index to end read variable from netcdf in additional dims
-    INTEGER(iwp)                                     :: z_start                         !> Index to start read variable from netcdf in additional dims
-    INTEGER(iwp)                                     :: z_end                           !> Index to end read variable from netcdf in additional     dims
-    INTEGER(iwp)                                     :: dt_emis                         !> Time Step Emissions 
-    INTEGER(iwp)                                     :: len_index                       !> length of index (used for several indices)
-    INTEGER(iwp)                                     :: len_index_voc                   !> length of voc index
-    INTEGER(iwp)                                     :: len_index_pm                    !> length of PMs index
-    REAL(wp)                                         :: con_factor                      !> Units Conversion Factor
- 
-
-    ! ---------------------------------------------------------------
-    ! Set Values of molecules, mols, etc
-    ! ---------------------------------------------------------------
-  
-    !> Avogadro number
-    REAL, PARAMETER        ::  Avog = 6.02205e23    ! mlc/mol
- 
-    !> Dobson units:
-    REAL, PARAMETER        ::  Dobs = 2.68668e16    ! (mlc/cm2) / DU
-
-    !> sesalt composition:
-    ! (Seinfeld and Pandis, "Atmospheric Chemistry and Physics",
-    !  table 7.8 "Composition of Sea-Salt", p. 444)
-    REAL, PARAMETER        ::  massfrac_Cl_in_seasalt  = 0.5504       ! (kg Cl )/(kg seasalt)
-    REAL, PARAMETER        ::  massfrac_Na_in_seasalt  = 0.3061       ! (kg Na )/(kg seasalt)
-    REAL, PARAMETER        ::  massfrac_SO4_in_seasalt = 0.0768       ! (kg SO4)/(kg seasalt)
-  
-    !> other numbers 
-    REAL, PARAMETER        ::  xm_seasalt =  74.947e-3                ! kg/mol : NaCl, SO4, ..
-    REAL, PARAMETER        ::  rho_seasalt = 2.2e3                    ! kg/m3
-
-    !> * amonium sulphate
-  
-    REAL, PARAMETER        ::  xm_NH4HSO4  =  xm_NH4 + xm_H + xm_SO4  ! kg/mol
-    REAL, PARAMETER        ::  rho_NH4HSO4a = 1.8e3                   ! kg/m3
-
-    ! ---------------------------------------------------------------
-    ! gas
-    ! ---------------------------------------------------------------
-  
-    !> gas constant                        
-    REAL, PARAMETER        ::  Rgas = 8.3144     ! J/mol/K
-  
-    !> gas constant for dry air
-    REAL, PARAMETER        ::  Rgas_air = Rgas / xm_air   ! J/kg/K
-  
-    !> water vapour
-    REAL, PARAMETER        ::  Rgas_h2o = Rgas / xm_h2o   ! J/kg/K
-  
-    !> standard pressure
-    REAL, PARAMETER        ::  p0 = 1.0e5    ! Pa
-
-    !> sea-level pressure
-    REAL, PARAMETER        ::  p_sealevel = 1.01325e5    ! Pa  <-- suggestion Bram Bregman
-
-    !> global mean pressure
-    REAL, PARAMETER        ::  p_global = 98500.0   ! Pa
-
-    !> specific heat of dry air at constant pressure
-    REAL, PARAMETER        ::  cp_air = 1004.0           ! J/kg/K
-
-    !> Latent heat of evaporation
-    REAL, PARAMETER        ::  lvap = 2.5e6     !  [J kg-1]
-
-    !> Latent heat of condensation at 0 deg Celcius
-    !  (heat (J) necesarry to evaporate 1 kg of water)
-    REAL, PARAMETER        ::  Lc = 22.6e5           ! J/kg
-  
-    !> kappa = R/cp = 2/7
-    REAL, PARAMETER        ::  kappa = 2.0/7.0
-
-    !> Von Karman constant (dry_dep)
-    REAL, PARAMETER        ::  vkarman = 0.4
-
-    !> Boltzmann constant:
-    REAL(wp), PARAMETER     ::  kbolz = 1.38066e-23_wp    ! J/K
-
-    !> Inverse Reference Pressure (1/Pa)    
-    REAL(wp), PARAMETER     ::  pref_i  = 1.0_wp / 100000.0_wp       
- 
-    !>
-    REAL(wp), PARAMETER     ::  r_cp    = 0.286_wp                  !< R / cp (exponent for potential temperature)
+    CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
+
+    INTEGER(iwp) ::  i                               !< index 1st selected dimension (some dims are not spatial)
+    INTEGER(iwp) ::  j                               !< index 2nd selected dimension 
+    INTEGER(iwp) ::  i_start                         !< Index to start read variable from netcdf along one dims
+    INTEGER(iwp) ::  i_end                           !< Index to end read variable from netcdf in one dims
+    INTEGER(iwp) ::  j_start                         !< Index to start read variable from netcdf in additional dims
+    INTEGER(iwp) ::  j_end                           !< Index to end read variable from netcdf in additional dims
+    INTEGER(iwp) ::  z_start                         !< Index to start read variable from netcdf in additional dims
+    INTEGER(iwp) ::  z_end                           !< Index to end read variable from netcdf in additional dims
+    INTEGER(iwp) ::  dt_emis                         !< Time Step Emissions 
+    INTEGER(iwp) ::  len_index                       !< length of index (used for several indices)
+    INTEGER(iwp) ::  len_index_voc                   !< length of voc index
+    INTEGER(iwp) ::  len_index_pm                    !< length of PMs index
+
+    REAL(wp) ::  con_factor                          !< Units Conversion Factor
+                           
+    REAL(wp), PARAMETER ::  Rgas = 8.3144                    !< gas constant in J/mol/K            
+    REAL(wp), PARAMETER ::  pref_i  = 1.0_wp / 100000.0_wp   !< Inverse Reference Pressure (1/Pa) 
+    REAL(wp), PARAMETER ::  r_cp    = 0.286_wp               !< R / cp (exponent for potential temperature)
 
 
     SAVE
 
-!-- Interfaces Part
-
+
+!    
 !-- Checks Input parameters
     INTERFACE chem_emissions_check_parameters
        MODULE PROCEDURE chem_emissions_check_parameters
     END INTERFACE chem_emissions_check_parameters
-
-!-- Reading of NAMELIST parameters
-!    INTERFACE chem_emissions_parin
-!       MODULE PROCEDURE chem_emissions_parin
-!    END INTERFACE chem_emissions_parin
-
-!-- Output of information to the header file
-!    INTERFACE chem_emissions_header
-!       MODULE PROCEDURE chem_emissions_header
-!    END INTERFACE chem_emissions_header
-
+!
 !-- Matching Emissions actions  
     INTERFACE chem_emissions_match
        MODULE PROCEDURE chem_emissions_match
     END INTERFACE chem_emissions_match
-
+!
 !-- Initialization actions  
     INTERFACE chem_emissions_init
        MODULE PROCEDURE chem_emissions_init
     END INTERFACE chem_emissions_init
-
+!
 !-- Setup of Emissions
     INTERFACE chem_emissions_setup
@@ -246,10 +178,9 @@
 
 
-!-- Public Interfaces
-    PUBLIC chem_emissions_init,chem_emissions_match,chem_emissions_setup
-
+    
+    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
+!
 !-- Public Variables
-
-    PUBLIC con_factor, len_index,len_index_voc,len_index_pm
+    PUBLIC con_factor, len_index, len_index_pm, len_index_voc
 
  CONTAINS
@@ -263,53 +194,22 @@
 
 
-    !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?
-
     IMPLICIT NONE
 
-    INTEGER(iwp)           ::  tmp
-
-    TYPE(chem_emis_att_type)   ::  emt
-!
-!--    Call routine for controlling chem_emissions namelist 
-!    CALL chem_emissions_parin
-
-!TBD: In case the given emission values do not coincide with the passed names we should think of a solution: 
-!  I would personally do that if the name of the species differ from the number of emission values, I would 
-!  print a warning that says that in correspondance of that particular species the value is zero. 
-!  An alternative would be to initialize the cs_surface_flux values to a negative number. 
-
-!-- Check Emission Species Number equal to number of passed names for the chemistry species: 
-   IF ( SIZE(emt%species_name) /= emt%nspec  )  THEN
-
-       message_string = 'Numbers of input emission species names and number of species'             //          &
+    INTEGER(iwp) ::  tmp
+
+    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
+
+       message_string = 'Numbers of input emission species names and number of species'     //      &
                          'for which emission values are given do not match'                  
        CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 ) 
-
-
+            
     ENDIF
-
-    !-- Check Emission Species Number equals to number of passed names for the species
-    !IF ( SIZE(emt%species_name) /= SIZE(emt%species_index)  )  THEN
-    !      message_string = 'Number of input emission species names and '             //          &
-    !                       ' number of passed species indices do not match'                  
-    !      CALL message( 'chem_emissions_check_parameters', 'CM0101', 2, 2, 0, 6, 0 ) 
-
-    !ENDIF
-
-    !-- Check Emission Categories
-!    IF ( SIZE(chem_emis%cat_name) /= SIZE(chem_emis%cat_index)  )  THEN
-!       WRITE( message_string, * )                                                        &
-!       'Number of input emission categories name and categories index does not match\\'  
-!       CALL message( 'chem_emissions_check_parameters', 'CM0101', 1, 2, 0, 6, 0 ) 
-!    ENDIF
-
-
-
-    !TBD: Check which other consistency controls should be included
-
-   !TBD: Include check for spatial dimension of netcdf file. If they do not match with the ones 
-   !     of the simulation, the model should print an error.  
-
-END SUBROUTINE chem_emissions_check_parameters
+    
+
+ END SUBROUTINE chem_emissions_check_parameters
 
 !------------------------------------------------------------------------------!
@@ -317,328 +217,352 @@
 ! ------------
 !> 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.
+!> 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
+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
-
-
-    !> Tell the user where we are!!
+    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
+
+
     CALL location_message( 'Matching input emissions and model chemistry species', .FALSE. )
 
-    !> Number of input emission species.
-
+    !
+    !-- Number of input emission species.
     nspec_emis_inp=emt_att%nspec 
 
-    !> Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED  !TBD: Add option for capital or small letters
-    SELECT CASE(TRIM(mode_emis)) 
-
-       !> PRE-PROCESSED case: in this case the input species have to coincide with the ones of the model (except VOCs for which we have the VOC split): NO and NO2 in input and not NOx.
-       CASE ("PRE-PROCESSED")
-
-          CALL location_message( 'chem_emissions PRE-PROCESSED_mode_matching', .FALSE. )
-
-          len_index=0
-          len_index_voc=0 
-
-          !> The first condition is that both the number of model and input emissions species are not null 
-          IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0)) THEN
-
-             !> Cycle over model species
-             DO ind_mod = 1, SIZE(spc_names)
-                !> Cycle over Input species  
+    !
+    !-- Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED 
+    SELECT CASE( TRIM( mode_emis ) ) 
+
+       !
+       !-- PRE-PROCESSED mode
+       CASE ( "PRE-PROCESSED" )
+
+          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
 
-           !> In the PRE-PROCESSED mode each emission species is given input values, made exception for the VOCs, having the total number of VOCs in input: the model then executes a split through the different VOC species
-                   ! > Check for VOC Species: In input in this case we only have one species (VOC)  
-                   IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN
-                      !> Cycle over the input voc species: they have to be one of the VOCs of the mechanism used. A list of VOC species for different mechanisms is provided in the module documentation
-                      DO ind_voc= 1,emt_att%nvoc
-                         !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays 
-                         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
+                   !
+                   !-- 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
+                   !
+                   !-- 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
-
-                   ALLOCATE(match_spec_voc_model(len_index_voc))   !> Contains indices of the VOC model species 
-
-                   ALLOCATE(match_spec_voc_input(len_index_voc))   !> In input there is only one value for VOCs in the DEFAULT mode. This array contains the indices of the different values of VOC compositions of the input variable VOC_composition
+             !
+             !-- 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
 
-                !> Repeat the same cycle of above but passing the species indices to the newly declared arrays
-                len_index=0
-
-                !> Cycle over model species
-                DO ind_mod = 1, SIZE(spc_names)  
-                   !> Cycle over Input species  
+                !
+                !-- 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
-
-                   !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays 
-
-                      ! > 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
+                      !
+                      !-- 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                         
+                               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
+                      !
+                      !-- 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
 
-             !> In the case there are no species matching, the emission module should not be called
              ELSE
-
-                message_string = 'Given Chemistry Emission Species'            //          &
-                                 ' DO NOT MATCH'                                //          &
-                                 ' model chemical species:'                      //          &
-                                 ' Chemistry Emissions routine is not called'  
+                !
+                !-- 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 !> IF (len_index>0)
+             ENDIF
  
           ELSE
 
-             !In this case either spc_names is null or nspec_emis_inp is not allocated
-
+             !
+             !-- 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:'                //          &
-                              ' Chemistry Emissions routine is not called'                  
+                              ' 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 ( (SIZE(spc_names) .GT. 0) .AND. ALLOCATED(nspec_emis_inp))
-
-       !> ------------------------------------------------------------------
-       !> DEFAULT case
-
+          ENDIF
+
+       !
+       !-- DEFAULT mode
        CASE ("DEFAULT")
 
-          CALL location_message( 'chem_emissions DEFAULT_mode_matching', .FALSE. )
-
-          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 PM Types: PM1, PM2.5, PM10. It is needed because the number of emission inputs and the one of PMs in the model may not be the same.
-
-          !> The first condition is that both the number of model and input emissions species are not null 
-          IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN
-
-             !> Cycle over model species
-             DO ind_mod = 1, SIZE(spc_names)
-                !> Cycle over input species
+          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: In input in this case we only have one species (VOC)  
-                   IF (TRIM(emt_att%species_name(ind_inp))=="VOC") THEN
-                      !> Cycle over the voc species given in input: they have to be one of the VOCs of the mechanism used. A list of VOC species for different mechanisms is provided in the module documentation 
-                      DO ind_voc= 1,emt_att%nvoc
-                         !> Determine the indices of the coinciding species in the two cases and assign them to matching arrays 
-                         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
+                   !
+                   !-- 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: For PMs there is only one input species name for all the PM types. This variable has 3 dimensions, one for each PM type.
-                   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
-                      ELSE IF (TRIM(spc_names(ind_mod))=="PM25") THEN
-                         len_index=len_index+1
-                      !>pm10
-                      ELSE IF (TRIM(spc_names(ind_mod))=="PM10") THEN
-                         len_index=len_index+1   
+                   !
+                   !-- 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: for NOx by definition we have 2 splits. The Emission Input Name in this case is only one: NOx, while in the model we can have 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
-                      ELSE IF (TRIM(spc_names(ind_mod))=="NO2") THEN
-                         len_index=len_index+1                        
+                   !
+                   !-- 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: same as for NOx, but with SO2 and SO4
-                   IF (TRIM(emt_att%species_name(ind_inp))=="SOx") THEN
-                      !>no
-                      IF (TRIM(spc_names(ind_mod))=="SO2") THEN
-                         len_index=len_index+1
-                      !>no2
-                      ELSE IF (TRIM(spc_names(ind_mod))=="SO4") THEN
-                         len_index=len_index+1                        
+                   !
+                   !-- 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
+                   !
+                   !-- Other Species
+                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                      len_index = len_index + 1
                    ENDIF
-                END DO   !> number of emission input species
-             END DO   !> number of model species
-
-
-             !> 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
-                   ALLOCATE(match_spec_voc_model(len_index_voc))   !> contains indices of the VOC model species 
-                   ALLOCATE(match_spec_voc_input(len_index_voc))   !> In input there is only one value for VOCs in the DEFAULT mode. 
-                                                                   !  This array contains the indices of the different values of VOC compositions of the input variable VOC_composition
+                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
 
-                !> ASSIGN INDICES
-                !> Repeat the same cycles of above, but passing the species indices to the newly declared arrays
-                len_index=0
-                len_index_voc=0
+                !
+                !-- Pass the species indices to declared arrays
+                len_index = 0
+                len_index_voc = 0
                 
-                DO ind_mod = 1, SIZE(spc_names)
+                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
+                      !
+                      !-- 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                         
+                               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:In this case the Inputs have only PM while the model has three different possible types: PM1, PM2.5, PM10. We need an additional index for matching each PM index in the model. 
-                      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
+                      !
+                      !-- 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
   
-                            !match_spec_pm(1)=ind_mod
-                         !>pm2.5
-                         ELSE IF (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
-  
-                            !match_spec_pm(2)=ind_mod
-                         !>pm10
-                         ELSE IF (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
-  
-                            !match_spec_pm(3)=ind_mod
                          ENDIF
                       ENDIF
 
-                      !> NOx - The same as for PMs but here the model species are only 2: NO and NO2
-                      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
+                      !
+                      !-- 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
   
-                            !match_spec_nox(1)=ind_mod
-                         !>no2
-                         ELSE IF (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
-  
-                           ! match_spec_nox(2)=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
+                      !
+                      !-- 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
   
-                           ! match_spec_sox(1)=ind_mod
-                         !>so4
-                         ELSE IF (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
+                         !
+                         !-- 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
   
-                           ! match_spec_sox(2)=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
+                      !
+                      !-- 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
+                         match_spec_input(len_index) = ind_inp
+                         match_spec_model(len_index) = ind_mod
                       ENDIF
                    END DO
@@ -647,8 +571,8 @@
              ELSE
 
-                message_string = 'Given Chemistry Emission Species'            //          &
-                                 ' DO NOT MATCH'                                //          &
-                                 ' model chemical species'                      //          &
-                                 ' Chemistry Emissions routine is not called'         
+                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 ) 
 
@@ -658,63 +582,63 @@
 
              message_string = 'Array of Emission species not allocated: '            //          &
-                              ' Either no emission species are provided as input or'  //          &
-                              ' no chemical species are used by PALM:'                //          &
-                              ' Chemistry Emissions routine is not called'                                   
+                              ' 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
-
-!-- CASE parameterized: In the parameterized case the user gives in input the exact species names of the chemical mechanism: no split is required for NOx, SOx, PMs and VOCs, but units of emissions inputs must be in (mole/m**)/s, made exception for PMs. 
-
+ 
+       !
+       !-- PARAMETERIZED mode
        CASE ("PARAMETERIZED")
 
-          len_index=0
-
-         !spc_names have to be greater than zero for proceeding
-          IF ( (SIZE(spc_names) .GT. 0) .AND. (nspec_emis_inp .GT. 0) ) THEN
-
-
-             !cycle over model species
-             DO ind_mod = 1, SIZE(spc_names)
-                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
-                      len_index=len_index+1
+          len_index = 0
+
+          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
+                      len_index = len_index + 1
                    ENDIF
-                   ind_inp=ind_inp+1
-
+                   ind_inp = ind_inp + 1
                 ENDDO
              ENDDO
 
-             !Proceed only if there are matched species
-             IF ( len_index .GT. 0 ) THEN
-
-                !Allocation of Arrays of the matched species
-                ALLOCATE(match_spec_input(len_index)) 
- 
-                ALLOCATE(match_spec_model(len_index))
-
-                !Assign corresponding indices of input and model matched species
-                len_index=0
-
-                DO ind_mod = 1, SIZE(spc_names)  
+             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    
+                len_index = 0
+
+                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
-                         len_index=len_index+1
-                         match_spec_input(len_index)=ind_inp
-                         match_spec_model(len_index)=ind_mod
+                   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
+                         match_spec_model(len_index) = ind_mod
                       ENDIF
-                   ind_inp=ind_inp+1
+                   ind_inp = ind_inp + 1
                    END DO
                 END DO
 
-                ! also, Add AN if to check that when the surface_csflux are passed to the namelist, also the street factor values are provided
-                DO  ispec = 1 , len_index
-
-                   IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. &
-                        emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
-
-                      message_string = 'PARAMETERIZED emissions mode selected:'            //          &
+                !
+                !-- Check
+                DO ispec = 1, len_index
+
+                   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:'            //           &
                                        ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
                                        ' but values of scaling factors for street types'    //          &
@@ -730,8 +654,8 @@
              ELSE
                 
-                message_string = 'Given Chemistry Emission Species'            //          &
-                                 ' DO NOT MATCH'                                //          &
-                                 ' model chemical species'                      //          &
-                                 ' Chemistry Emissions routine is not called'         
+                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 ) 
              ENDIF
@@ -739,8 +663,8 @@
           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.'                //          &
-                              ' Chemistry Emissions routine is not called'                                   
+                              ' Emission routine is not called'                                   
              CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) 
  
@@ -748,8 +672,9 @@
 
 
-!-- CASE when emission module is switched on but mode_emis is not specified or it is given the wrong name
+       !
+       !-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
        CASE DEFAULT
 
-          message_string = 'Chemistry Emissions Module switched ON, but'            //          &
+          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"'                  
@@ -758,118 +683,117 @@
        END SELECT
 
-END SUBROUTINE chem_emissions_match
-
+ END SUBROUTINE chem_emissions_match
+
+ 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Initialization:
-!> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0
-!>  
+!> Netcdf reading, arrays allocation and first calculation of cssws
+!> fluxes at timestep 0 
 !------------------------------------------------------------------------------!
- SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out)
+ SUBROUTINE chem_emissions_init( emt_att, emt, nspec_out )
 
     USE surface_mod,                                                           &
-       ONLY:  surf_lsm_h,surf_def_h,surf_usm_h
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
     
     IMPLICIT NONE
  
-    TYPE(chem_emis_att_type),INTENT(INOUT)                            :: emt_att   !> variable where to store all the information of 
-                                                                                   !  emission inputs whose values do not depend
-                                                                                   !  on the considered species
-
-    TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt      !> variable where to store emission inputs values, 
-                                                                                   !  depending on the considered species
+    CHARACTER (LEN=80)          ::  units                                           !< units of inputs
     
-    INTEGER(iwp),INTENT(INOUT)                                        :: nspec_out !> number of outputs of chemistry emission routines
-
-    CHARACTER (LEN=80)                                                :: units     !> units of chemistry inputs
- 
-    INTEGER(iwp)                                                      :: ispec     !> Index to go through the emission chemical species
-
-
-!-- Actions for initial runs : TBD: needs to be updated
-  IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+    INTEGER(iwp), INTENT(INOUT) ::  nspec_out                                       !< number of outputs
+
+    INTEGER(iwp)                :: ispec                                            !< running index
+
+    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att                             !< variable where to store all the information
+                                                                                    !< of emission inputs whose values do not 
+                                                                                    !< depend on the considered species
+
+    TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt      !< variable where to store emission input  
+                                                                                    !< values depending on the considered species
+
+!    
+!-- Actions for initial runs 
+!  IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 !--    ... 
-    
+!    
 !
 !-- Actions for restart runs
+!  ELSE
+!--    ...
+!
+!  ENDIF
+
+
+  CALL location_message( 'Starting initialization of emission arrays', .FALSE. )
+
+  !
+  !-- Matching
+  CALL  chem_emissions_match( emt_att, nspec_out )
+ 
+  IF ( nspec_out == 0 )  THEN
+  
+     emission_output_required = .FALSE.
+
   ELSE
-!--    ...
-
-  ENDIF
-
-
-  CALL location_message( 'Starting initialization of chemistry emissions arrays', .FALSE. )
-
-
-  !-- Match Input and Model emissions
-  CALL  chem_emissions_match(emt_att,nspec_out)
-
-  !> If nspec_out==0, then do not use emission module: The corresponding message is already printed in the matching routine 
-  IF ( nspec_out == 0 ) THEN
-  
-     emission_output_required=.FALSE.
-
-  ELSE
-
-     emission_output_required=.TRUE.
-
-
-!-----------------------------------------------------------------
-     !> Set molecule masses' 
-     ALLOCATE(emt_att%xm(nspec_out))
-     ! loop over emisisons:
+
+     emission_output_required = .TRUE.
+
+
+     !
+     !-- Set molecule masses' 
+     ALLOCATE( emt_att%xm(nspec_out) )
 
      DO ispec = 1, nspec_out
-       ! switch:
-        SELECT CASE ( TRIM(spc_names(match_spec_model(ispec))) )
-           CASE ( 'SO2','SO4' ) ; emt_att%xm(ispec) = xm_S + xm_O * 2      ! kg/mole
-           CASE ( 'NO','NO2' )  ; emt_att%xm(ispec) = xm_N + xm_O * 2      ! kg/mole  NO2 equivalent
-           CASE ( 'NH3' ) ; emt_att%xm(ispec) = xm_N + xm_H * 3  ! kg/mole
-           CASE ( 'CO'  ) ; emt_att%xm(ispec) = xm_C + xm_O      ! kg/mole
-           CASE ( 'CO2' ) ; emt_att%xm(ispec) = xm_C + xm_O * 2  ! kg/mole
-           CASE ( 'CH4' ) ; emt_att%xm(ispec) = xm_C + xm_H * 4  ! kg/mole     
-           CASE ( 'HNO3' ); emt_att%xm(ispec) = xm_H + xm_N + xm_O*3 !kg/mole  
+        SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
+           CASE ( 'SO2' ); emt_att%xm(ispec) = xm_S + xm_O * 2        !< kg/mole
+           CASE ( 'SO4' ); emt_att%xm(ispec) = xm_S + xm_O * 4        !< kg/mole
+           CASE ( 'NO' ); emt_att%xm(ispec) = xm_N + xm_O             !< kg/mole
+           CASE ( 'NO2' ); emt_att%xm(ispec) = xm_N + xm_O * 2        !< kg/mole   
+           CASE ( 'NH3' ); emt_att%xm(ispec) = xm_N + xm_H * 3        !< kg/mole
+           CASE ( 'CO'  ); emt_att%xm(ispec) = xm_C + xm_O            !< kg/mole
+           CASE ( 'CO2' ); emt_att%xm(ispec) = xm_C + xm_O * 2        !< kg/mole
+           CASE ( 'CH4' ); emt_att%xm(ispec) = xm_C + xm_H * 4        !< kg/mole     
+           CASE ( 'HNO3' ); emt_att%xm(ispec) = xm_H + xm_N + xm_O*3  !< kg/mole  
            CASE DEFAULT
-             !--  TBD: check if this hase to be removed
-              !emt_att%xm(ispec) = 1.0_wp
+              emt_att%xm(ispec) = 1.0_wp
         END SELECT
      ENDDO
 
     
-     !> ASSIGN EMISSION VALUES for the different emission modes
-     SELECT CASE(TRIM(mode_emis))   !TBD: Add the option for CApital or small letters 
-
-
-        !> PRE-PROCESSED case
-        CASE ("PRE-PROCESSED")
-
-           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out))  
-
-           CALL location_message( 'emis_distribution array allocated in PRE-PROCESSED mode', .FALSE. )
- 
-           !> Calculate the values of the emissions at the first time step
-           CALL chem_emissions_setup(emt_att,emt,nspec_out)
-
-        !> Default case
-        CASE ("DEFAULT")
-
-           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution( 1, 0:ny, 0:nx, nspec_out ) ) 
-        
-           CALL location_message( 'emis_distribution array allocated in DEFAULT mode', .FALSE. )
-
-           !> Calculate the values of the emissions at the first time step
-           CALL chem_emissions_setup(emt_att,emt,nspec_out)
-
-        !> PARAMETERIZED case
-        CASE ("PARAMETERIZED")
-
-           CALL location_message( 'emis_distribution array allocated in PARAMETERIZED mode', .FALSE. )
-
-           ! For now for PAR and DEF values only, first vertical level of emis_distribution is allocated, while for PRE-PROCESSED all. 
-           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE(emis_distribution(1,0:ny,0:nx,nspec_out))
-
-           !> Calculate the values of the emissions at the first time step
-           CALL chem_emissions_setup(emt_att,emt,nspec_out)
+     !
+     !-- 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( emt_att, emt, 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( emt_att, emt, 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( emt_att, emt, nspec_out)
 
      END SELECT
@@ -887,99 +811,78 @@
 !-------------------------------------------------------------------------------!
 
- SUBROUTINE chem_emissions_setup(emt_att,emt,nspec_out)
+ SUBROUTINE chem_emissions_setup( emt_att, emt, nspec_out )
 
    USE arrays_3d,                                                    &
        ONLY:  dzw
-   USE grid_variables,                                                        &
+   USE grid_variables,                                               &
        ONLY: dx, dy
-   USE indices,                                                               &
-       ONLY: nnx,nny,nnz 
-   USE surface_mod,                                                           &
-       ONLY:  surf_lsm_h,surf_def_h,surf_usm_h
-   USE netcdf_data_input_mod,                                                 &
+   USE indices,                                                      &
+       ONLY: nnx, nny, nnz 
+   USE surface_mod,                                                  &
+       ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+   USE netcdf_data_input_mod,                                        &
        ONLY:  street_type_f
-   USE arrays_3d,                                                             &        
+   USE arrays_3d,                                                    &        
        ONLY: hyp, pt 
 
+   
  IMPLICIT NONE
   
-    !--- IN/OUT
-
-    TYPE(chem_emis_att_type),INTENT(INOUT)                            ::  emt_att    !> variable where to store all the information of 
-                                                                                     !  emission inputs whose values do not depend
-                                                                                     !  on the considered species
-
-    TYPE(chem_emis_val_type),INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt        !> variable where to store emission inputs values, 
-                                                                                     !  depending on the considered species
-
-    INTEGER,INTENT(IN)                                                ::  nspec_out  !> Output of routine chem_emis_matching with number
-                                                                                     !  of matched species
-    !---
-
-    REAL(wp),ALLOCATABLE,DIMENSION(:,:)                               ::  delta_emis !> Term to add to the emission_distribution array 
-                                                                                     !  for each of the categories at each time step 
-                                                                                     !  in the default mode
-
-    CHARACTER(LEN=80)                                                 ::  units      !> Units of the emissions
-
-    INTEGER(iwp)                                                      ::  icat       !> Index for number of categories
-    INTEGER(iwp)                                                      ::  ispec      !> index for number of species
-    INTEGER(iwp)                                                      ::  i_pm_comp  !> index for number of PM components
-    INTEGER(iwp)                                                      ::  ivoc       !> Index for number of components of  VOCs
-    INTEGER(iwp)                                                      ::  time_index !> Index for time
+
+    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
+                                                                                !< of matched species
+
+    CHARACTER(LEN=80) ::  units                                                 !< Units of the emission data
+
+    INTEGER(iwp) ::  i                                                          !< running index for grid in x-direction
+    INTEGER(iwp) ::  j                                                          !< running index for grid in y-direction
+    INTEGER(iwp) ::  k                                                          !< running index for grid in z-direction
+    INTEGER(iwp) ::  m                                                          !< running index for horizontal surfaces
     
-    REAL(wp),ALLOCATABLE, DIMENSION(:)                                ::  time_factor !> factor for time scaling of emissions
-    REAL(wp),ALLOCATABLE, DIMENSION(:,:)                              ::  emis
-
-    REAL(wp), DIMENSION(24)                                           :: par_emis_time_factor      !< time factors 
-                                                                                      !  for the parameterized mode: these are fixed for each hour
-                                                                                      !  of a single day.
-    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)                    ::  conv_to_ratio !> factor used for converting input 
-                                                                                        !  to adimensional concentration ratio
-
-   ! ------------------------------------------
-    ! variables for the conversion of indices i and j according to surface_mod
-
-    INTEGER(iwp)                                                      ::  i          !> running index for grid in x-direction
-    INTEGER(iwp)                                                      ::  j          !> running index for grid in y-direction
-    INTEGER(iwp)                                                      ::  k          !> running index for grid in z-direction
-    INTEGER(iwp)                                                      ::  m          !> running index for horizontal surfaces
-
-    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)                    ::  tmp_temp
-
-    ! --- const -------------------------------
-    !-CONVERSION FACTORS: TIME
-    ! number of sec per hour = 3600   
-    REAL, PARAMETER   ::  s_per_hour = 3600.0  !  (s)/(hour)
-    ! number of sec per day = 86400   
-    REAL, PARAMETER   ::  s_per_day = 86400.0  !  (s)/(day)
-    ! number of hours in a year of 365 days:
-    REAL, PARAMETER   ::  hour_per_year = 8760.0 !> TBD: What about leap years?
-    ! number of hours in a day:
-    REAL, PARAMETER   ::  hour_per_day = 24.0
-
-    ! conversion from hours to seconds (s/hour) = 1/3600.0 ~ 0.2777778   
-    REAL, PARAMETER   ::  hour_to_s = 1/s_per_hour  !  (hour)/(s)
-    ! conversion from day to seconds (s/day) = 1/86400.0 ~ 1.157407e-05   
-    REAL, PARAMETER   ::  day_to_s = 1/s_per_day   !  (day)/(s)
-    ! conversion from year to sec (s/year) = 1/31536000.0 ~ 3.170979e-08   
-    REAL, PARAMETER   ::  year_to_s = 1/(s_per_hour*hour_per_year)  !  (year)/(s)
-
-    !-CONVERSION FACTORS: WEIGHT
-    !  Conversion from tons to kg (kg/tons) = 100.0/1 ~ 100
-    REAL, PARAMETER   ::  tons_to_kg = 100  !  (tons)/(kg)
-    !  Conversion from g to kg (kg/g) = 1/1000.0 ~ 0.001
-    REAL, PARAMETER   ::  g_to_kg = 0.001  !  (g)/(kg)
-    !  Conversion from g to kg (kg/g) = 1/1000.0 ~ 0.001
-    REAL, PARAMETER   ::  miug_to_kg = 0.000000001  !  (microng)/(kg)
-
-    !< CONVERSION FACTORS: fraction to ppm
-    REAL(wp), PARAMETER   ::  ratio2ppm  = 1.0e06_wp  
+    INTEGER(iwp) ::  icat                                                       !< Index for number of categories
+    INTEGER(iwp) ::  ispec                                                      !< index for number of species
+    INTEGER(iwp) ::  i_pm_comp                                                  !< index for number of PM components
+    INTEGER(iwp) ::  ivoc                                                       !< Index for number of VOCs
+    INTEGER(iwp) ::  time_index                                                 !< Index for time
+
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  delta_emis                        
+    REAL(wp), ALLOCATABLE, DIMENSION(:)   ::  time_factor                       !< factor for time scaling of emissions
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  emis
+
+    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 
+                                                                                !< to concentration ratio
+    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  tmp_temp
+
+    !
+    !-- CONVERSION FACTORS: TIME  
+    REAL(wp), PARAMETER ::  s_per_hour = 3600.0                       !< number of sec per hour (s)/(hour)    
+    REAL(wp), PARAMETER ::  s_per_day = 86400.0                       !< number of sec per day (s)/(day)  
+    REAL(wp), PARAMETER ::  hour_per_year = 8760.0                    !< number of hours in a year of 365 days  
+    REAL(wp), PARAMETER ::  hour_per_day = 24.0                       !< number of hours in a day
+    
+    REAL(wp), PARAMETER ::  hour_to_s = 1/s_per_hour                  !< conversion from hours to seconds (s/hour) ~ 0.2777778      
+    REAL(wp), PARAMETER ::  day_to_s = 1/s_per_day                    !< conversion from day to seconds (s/day) ~ 1.157407e-05
+    REAL(wp), PARAMETER ::  year_to_s = 1/(s_per_hour*hour_per_year)  !< conversion from year to sec (s/year) ~ 3.170979e-08
+    !
+    !-- CONVERSION FACTORS: WEIGHT  
+    REAL(wp), PARAMETER ::  tons_to_kg = 100                          !< Conversion from tons to kg (kg/tons)    
+    REAL(wp), PARAMETER ::  g_to_kg = 0.001                           !< Conversion from g to kg (kg/g) 
+    REAL(wp), PARAMETER ::  miug_to_kg = 0.000000001                  !< Conversion from g to kg (kg/g)
+    !
+    !-- CONVERSION FACTORS: fraction to ppm
+    REAL(wp), PARAMETER ::  ratio2ppm  = 1.0e06  
     !------------------------------------------------------    
 
-    IF ( emission_output_required ) THEN
-
-    !>  Set emis_dt  !TBD: this is the same as dt_chem. We should consider the fact that dt_emis should be the timestep of input emissions or better defined, the timestep at which the emission routines are called: for now one hour. It should be made changeable.
-
+    IF ( emission_output_required )  THEN
+
+       !
+       !-- Set emis_dt  
        IF ( call_chem_at_all_substeps )  THEN
 
@@ -993,88 +896,75 @@
 
 
-    ! --- CHECK UNITS
-    !>-----------------------------------------------------
-    !> Conversion of the units to the ones employed in PALM. 
-    !> Possible temporal units of input emissions data are: year, hour and second; 
-    !> the mass could be expressed in: tons, kilograms or grams. 
-    !> IN the PARAMETERIZED mode no conversion is possible: in this case INPUTS have to have fixed units.
-    !>-----------------------------------------------------
-
-        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') 
-                CALL location_message( 'Units of the emissions are consistent with'  //          &
-                                       ' the ones employed in the PALM-4U model ', .FALSE. )
+       !
+       !-- 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')
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
+             CASE ('kg/m2/hour', 'KG/M2/HOUR' )
 
                 con_factor=hour_to_s
 
-             CASE ( 'kg/m2/day','KG/M2/DAY' ) 
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'kg/m2/day', 'KG/M2/DAY' ) 
+                
                 con_factor=day_to_s
 
-             CASE ( 'kg/m2/year','KG/M2/YEAR' )
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'kg/m2/year', 'KG/M2/YEAR' )
+              
                 con_factor=year_to_s
 
-          !> Tons
-             CASE ( 'ton/m2/s','TON/M2/S' ) 
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             !
+             !-- Tons
+             CASE ( 'ton/m2/s', 'TON/M2/S' ) 
+               
                 con_factor=tons_to_kg
 
-             CASE ( 'ton/m2/hour','TON/M2/HOUR' ) 
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ) 
+               
                 con_factor=tons_to_kg*hour_to_s
 
-             CASE ( 'ton/m2/year','TON/M2/YEAR' ) 
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'ton/m2/year', 'TON/M2/YEAR' ) 
+               
                 con_factor=tons_to_kg*year_to_s
 
-          !> Grams
-             CASE ( 'g/m2/s','G/M2/S' )
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             !
+             !-- Grams
+             CASE ( 'g/m2/s', 'G/M2/S' )
+          
                 con_factor=g_to_kg
 
-             CASE ( 'g/m2/hour','G/M2/HOUR' ) 
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
+             CASE ( 'g/m2/hour', 'G/M2/HOUR' ) 
 
                 con_factor=g_to_kg*hour_to_s
 
-             CASE ( 'g/m2/year','G/M2/YEAR' )
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'g/m2/year', 'G/M2/YEAR' )
+ 
                 con_factor=g_to_kg*year_to_s
 
-          !> Micro Grams
-             CASE ( 'micrograms/m2/s','MICROGRAMS/M2/S' )
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             !
+             !-- Micrograms
+             CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' )
+ 
                 con_factor=miug_to_kg
 
-             CASE ( 'micrograms/m2/hour','MICROGRAMS/M2/HOUR' ) 
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ) 
+ 
                 con_factor=miug_to_kg*hour_to_s
 
-             CASE ( 'micrograms/m2/year','MICROGRAMS/M2/YEAR' )
-                CALL location_message( 'Units of emission inputs need conversion', .FALSE. )
-
+             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' )
+  
                 con_factor=miug_to_kg*year_to_s
 
              CASE DEFAULT  
-                message_string = 'The Units of the provided input chemistry emission species' // &
+                message_string = 'The Units of the provided emission input' // &
                                  ' are not the ones required by PALM-4U: please check '      // &
-                                 ' chemistry emission module documentation.'                                 
+                                 ' emission module documentation.'                                 
                 CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
 
@@ -1082,95 +972,88 @@
 
          
-       !> PRE-PROCESSED and parameterized mode
-       ELSE 
- 
-          message_string = 'No Units conversion required for units of chemistry emissions' // &
-                           ' of the PARAMETERIZED mode: units have to be provided in'     //  & 
-                           ' micromole/m**2/day for GASES and'                            //  &
-                           ' kg/m**2/day for PMs'                      
-          CALL message( 'chem_emissions_setup', 'CM0447', 0, 0, 0, 6, 0 ) 
-
        ENDIF
 
-    !> Conversion factors (if the units are kg/m**2/s) we have to convert them to ppm/s
-     DO i=nxl,nxr
-          DO j=nys,nyn
-          !> Derive Temperature from Potential Temperature
+       !
+       !-- Conversion factor to convert  kg/m**2/s to ppm/s
+       DO i = nxl, nxr
+          DO j = nys, nyn
+             !
+             !-- 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 basically we need to 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
 
-    !>------------------------------------------------
-    !> Start The Processing of the input data
-
-        emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
-
-    !>-----------------------------------------------
-    !> Distinguish between DEFAULT, PRE-PROCESSED and PARAMETERIZED mode when calculating time_factors: only done for DEFAULT mode. For the PARAMETERIZED mode there is a time factor but it is fixed in the model 
+
+       !
+       !-- Initialize
+       emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
+
   
-    !> PRE-PROCESSED MODE 
-       IF (TRIM(mode_emis)=="PRE-PROCESSED") THEN
-
-          !> Update time indices
-          CALL time_preprocessed_indices(index_hh)
-
-          CALL location_message( 'PRE-PROCESSED MODE: No time-factor specification required', .FALSE. )
-
-       ELSEIF (TRIM(mode_emis)=="DEFAULT") THEN
-
-          CALL location_message( 'DEFAULT MODE: time-factor specification required', .FALSE. )
- 
-       !> Allocate array where to store temporary emission values      
-          IF(.NOT. ALLOCATED(emis)) ALLOCATE(emis(nys:nyn,nxl:nxr))
-
-       !> Allocate time factor per emitted component 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 .LE. 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)      
+       !
+       !-- PRE-PROCESSED MODE 
+       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
+             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              
+                time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh)      
 
              ELSE
 
-                message_string = 'The "HOUR" time-factors (DEFAULT mode) of the chemistry emission species'           // &
+                message_string = 'The "HOUR" time-factors in the DEFAULT mode '           // &
                               ' are not provided for each hour of the total simulation time'      
                 CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 ) 
 
              ENDIF
-
-       !> 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) .LE. 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)*   &
+          !
+          !-- 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)
      
              ELSE
 
-                message_string = 'The "MDH" time-factors (DEFAULT mode) of the chemistry emission species'           // &
+                message_string = 'The "MDH" time-factors in the DEFAULT mode '           //  &
                               ' are not provided for each hour/day/month  of the total simulation time'      
                 CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
@@ -1179,27 +1062,26 @@
 
           ELSE
-          !> condition when someone used the DEFAULT mode but forgets to indicate the time-factor in the namelist
-             
-             message_string = 'The time-factor (DEFAULT mode) of the chemistry emission species'           // &
-                              ' is not provided in the NAMELIST'      
+                 
+             message_string = 'In the DEFAULT mode the time factor'           //  &
+                              ' has to be defined in the NAMELIST'      
              CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
          
           ENDIF
 
-       !> PARAMETERIZED MODE
-       ELSEIF (TRIM(mode_emis)=="PARAMETERIZED") THEN
-          CALL location_message( 'PARAMETERIZED MODE: time-factor specification is fixed: '                         // &
-                                 ' 24 values for every day of the year ', .FALSE. )
+       !
+       !-- PARAMETERIZED MODE
+       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, &
+          !
+          !-- 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 /)
           
-          !> in this case allocate time factor each hour in a day
-          IF (.NOT. ALLOCATED(time_factor)) ALLOCATE(time_factor(1))
-
-          !>Pass the values of time factors in the parameterized mode to the time_factor variable. in this case index_hh==hour_of_day
-          index_hh=hour_of_day
+          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)
@@ -1207,122 +1089,137 @@
        ENDIF
 
-!--------------------------------
-!--  EMISSION DISTRIBUTION Calculation
-
-       !> PARAMETERIZED CASE
-       IF ( mode_emis == "PARAMETERIZED" ) THEN
-
-          DO ispec=1,nspec_out
-
-             !> Values are still micromoles/(m**2*s). Units in this case are always 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
+
+       !
+       !-- 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
 
           ENDDO 
 
-       !> PRE-PROCESSED CASE
-       ELSEIF (TRIM(mode_emis)=="PRE-PROCESSED") THEN
-
-          !> Start Cycle over Species
-          DO ispec=1,nspec_out !> nspec_out represents the number of species in common between 
-                               !  the emission input data and the chemistry mechanism used
+       !
+       !-- 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
+             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
-
-      !TBD: At the moment the default mode considers a single vertical level (the surface). So we need to change it accordingly or eventually include the variable vertical_profile to be used in the default mode i we want to consider additional levels.
-
-       ELSE IF (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))  
-
-
-       !> Start Cycle over categories
+          
+       !
+       !-- 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
           DO icat = 1, emt_att%ncat
  
-          !> Start Cycle over Species
-             DO ispec=1,nspec_out !> nspec_out represents the number of species in common between 
-                                  !  the emission input data and the chemistry mechanism used
+             !
+             !-- 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)
 
-!TBD: The consideration of dt_emis of the input data is still missing. Basically the emissions could be provided every 10, 30 minutes and not always at one hour. This should be eventually solved in the date_and_time mode routine.
-
-             !> the time factors are 24 for each day. When multiplied by a daily value, they allow to have an hourly value. Then to convert it to seconds, we still have to divide this value by 3600. 
-             !  So given any units, we convert them to seconds and finally multiply them by 24 ((value/sec)*(24*3600)=value/day ---- (value/day)*time_factor=value/hour ---(value/hour)/(3600)=value/sec )
-             !                                                                                 ((value/sec)*(24*3600)*time_factor)/3600=24*(value/sec)*time_factor                          
-
-             !> NOX Compositions
-                IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN
-                !>             Kg/m2*s                   kg/m2*s
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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)
-
-                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN
-
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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)
- 
-             !> SOX Compositions
-
-                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN
-                   !>             Kg/m2*s                   kg/m2*s
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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)
-
-                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN
-                   !>             Kg/m2*s                   kg/m2*s
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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)
- 
-
-             !> PMs should be in kg/m**2/s, so conversion factor is here still required
-             !> PM1 Compositions
-                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1") THEN
-
-                !> Cycle over the different pm components for PM1 type
-                   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)* &
-                                                    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)+ &
+
+                !
+                !-- NOx
+                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) * 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) * 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)
                    ENDDO
 
-             !> PM2.5 Compositions
-                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM25") THEN
-
-                !> Cycle over the different pm components for PM2.5 type
-                   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)* &
-                                                    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
+
+                   !
+                   !-- 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)
  
                    ENDDO
 
-             !> PM10 Compositions
-                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
-
-                !> Cycle over the different pm components for PM10 type
-                   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)* &
-                                                    emt_att%pm_comp(icat,i_pm_comp,3)*con_factor*hour_per_day 
+                !
+                !-- 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)+ &
@@ -1331,17 +1228,18 @@
                    ENDDO
 
-             !> VOCs Compositions: for VOCs, the input value is provided in kg/m**2*s but the composition is provided in mole/kg: a conversion factor for the input that could be eventually provided in, for example, tons/(m**2*s) is still required
-                ELSE IF  (SIZE(match_spec_voc_input) .GT. 0) THEN
-
-                  !TBD: maybe we can avoid the cycle 
-                   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 )
+
+                      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                       
@@ -1349,39 +1247,40 @@
                    ENDDO
                 
-             !> General case (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  ! IF (spc_names==)
-
-                !> for every species and category set emis to 0 so to avoid overwriting. TBD: discuss whether necessary
-
+                   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
+                
                 emis(:,:)= 0 
-
+                
              ENDDO
-
-          !> for every category set delta_emis to 0 so to avoid overwriting. TBD: discuss whether necessary
-
-          delta_emis(:,:)=0 
- 
+             
+             delta_emis(:,:)=0 
+          
           ENDDO
 
-       ENDIF  !> mode_emis
-
-!-------------------------------------------------------------------------------------------------------
-!--- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
-!-------------------------------------------------------------------------------------------------------
-
+       ENDIF
+
+       
+!
+!-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
+!
 !-- PARAMETERIZED mode
-       !> For the PARAMETERIZED mode units of inputs are always micromoles/(m**2*s). In this case we do not need the molar mass for conversion into ppms
-       IF (TRIM(mode_emis)=="PARAMETERIZED") THEN
+       !
+       !-- 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
-             IF (surf_lsm_h%ns .GT. 0 ) THEN
+             !
+             !-- Streets are lsm surfaces, hence, no usm surface treatment required
+             IF ( surf_lsm_h%ns > 0 )  THEN
                 DO  m = 1, surf_lsm_h%ns
                    i = surf_lsm_h%i(m)
@@ -1389,34 +1288,33 @@
                    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 already matched species
-                      DO  ispec=1,nspec_out
-
-                         !> PMs are already in mass units: 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
+                      !
+                      !-- 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)) * &
-                            !                                                       kg/(m^2*s)
-                                                                                emis_distribution(1,j,i,ispec)*        &
-                            !                                                    kg/m^3
-                                                                                rho_air(k)   
-
+                                                                           emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
+                                                                            rho_air(k)                                    !< in kg/m^3
+                            
+                         !
+                         !-- Other Species
+                         !-- Inputs are micromoles
                          ELSE
 
-                         !> Other Species: inputs are micromoles: they have to be converted in moles
-                         !                 ppm/s *m *kg/m^3                
-                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*   &
-                         !                                                    micromoles/(m^2*s)
-                                                                          emis_distribution(1,j,i,ispec) *              &
-                         !                                                    m^3/Nmole
-                                                                          conv_to_ratio(k,j,i)*                         &       
-                         !                                                    kg/m^3
-                                                                          rho_air(k)   
- 
+                            !   
+                            !--             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
 
                          ENDIF
@@ -1427,30 +1325,30 @@
                    ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND. street_type_f%var(j,i) < main_street_id )  THEN
 
-                   !> Cycle over already matched species
-                      DO  ispec=1,nspec_out
-
-                         !> PMs are already in mass units: micrograms
-                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" &
-                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
-                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
-
-                            !              kg/(m^2*s) *kg/m^3                               
-                            surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *   &
-                            !                                                       kg/(m^2*s)
-                                                                                emis_distribution(1,j,i,ispec)*        &
-                            !                                                    kg/m^3
-                                                                                rho_air(k)   
+                      !
+                      !-- 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
                          ELSE
                        
-
-                         !>Other Species: inputs are micromoles
-                         !                 ppm/s *m *kg/m^3                
-                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) *   &
-                         !                                                    micromoles/(m^2*s)
-                                                                          emis_distribution(1,j,i,ispec) *              &
-                         !                                                    m^3/Nmole
-                                                                          conv_to_ratio(k,j,i)*                         &       
-                         !                                                    kg/m^3
-                                                                          rho_air(k)   
+                            !   
+                            !--             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   
                          ENDIF
 
@@ -1459,5 +1357,6 @@
                    ELSE
 
-                   !> If no street type is defined, then assign null emissions to all the species
+                      !
+                      !-- If no street type is defined, then assign zero emission to all the species
                       surf_lsm_h%cssws(:,m) = 0.0_wp
 
@@ -1465,16 +1364,19 @@
 
                 ENDDO
+
              ENDIF
 
           ENDIF
 
-       !> For both DEFAULT and PRE-PROCESSED
+       !
+       !-- For both DEFAULT and PRE-PROCESSED mode
        ELSE    
        
 
-          DO ispec=1,nspec_out 
-                                !TBD: for the PRE-PROCESSED mode in the future, values at different heights should be included!            
-             !> Default surface type
-             IF (surf_def_h(0)%ns .GT. 0) THEN
+          DO ispec = 1, nspec_out 
+                   
+             !
+             !-- Default surfaces
+             IF ( surf_def_h(0)%ns > 0 )  THEN
 
                 DO  m = 1, surf_def_h(0)%ns
@@ -1485,42 +1387,40 @@
                    IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
 
-
-                      !> Distinguish between PMs (no needing conversion in ppms), 
-                      !  VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
-                      ! other Species (still expressed in Kg/(m**2*s) at this point)
-
-                      !> 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
+                      !
+                      !-- 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)*   &
-                         !                                                  kg/m^3
-                                                                          rho_air(nzb)   
+                         !
+                         !--            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  
  
  
                       ELSE
 
-                         !> VOCs
-                         IF ( len_index_voc .GT. 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) *      &
-                                                                            !    m^3/mole          ppm
-                                                                             conv_to_ratio(nzb,j,i)*ratio2ppm *      &
-                         !                                                    kg/m^3
-                                                                             rho_air(nzb)   
-
-
-                         !> OTHER SPECIES
+                         !
+                         !-- 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
                          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) *      &
-                                                                               !  mole/Kg        
-                                                                             (1/emt_att%xm(ispec))*                &
-                            !                                                    m^3/mole          ppm
-                                                                             conv_to_ratio(nzb,j,i)*ratio2ppm*       &
-                            !                                                  kg/m^3
-                                                                             rho_air(nzb)   
+                            !
+                            !--           (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  
   
 
@@ -1533,8 +1433,9 @@
                 ENDDO
 
-             END IF
+             ENDIF
          
-             !-- Land Surface Mode surface type
-             IF (surf_lsm_h%ns .GT. 0) THEN
+             !
+             !-- LSM surfaces
+             IF ( surf_lsm_h%ns > 0 )  THEN
 
                 DO m = 1, surf_lsm_h%ns
@@ -1546,39 +1447,37 @@
                    IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN
 
-                      !> Distinguish between PMs (no needing conversion in ppms), 
-                      !  VOC (already converted to moles/(m**2*s) using conv_factors: they do not need molar masses for their conversion to PPMs ) and
-                      ! other Species (still expressed in Kg/(m**2*s) at this point)
-
-                      !> 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
+                      !
+                      !-- 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) *              &
-                         !                                                  kg/m^3
-                                                                          rho_air(k)   
+                         !
+                         !--         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
  
                       ELSE
 
-                         !> VOCs
-                         IF ( len_index_voc .GT. 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) *      &
-                                                                          !    m^3/mole          ppm
-                                                                          conv_to_ratio(k,j,i)*ratio2ppm*    &
-                         !                                                 kg/m^3
-                                                                          rho_air(k)   
-
-                         !> OTHER SPECIES
+                         !
+                         !-- 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
                          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) *               &
-                                                                               !  mole/Kg    
-                                                                         (1/emt_att%xm(ispec))*                          &
-                            !                                                m^3/mole           ppm
-                                                                         conv_to_ratio(k,j,i)*ratio2ppm*                 &
-                            !                                            kg/m^3
-                                                                         rho_air(k)   
+                            !
+                            !--           (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     
                                                       
                          ENDIF
@@ -1590,8 +1489,9 @@
                 ENDDO
 
-             END IF
-
-             !-- Urban Surface Mode surface type
-             IF (surf_usm_h%ns .GT. 0) THEN
+             ENDIF
+
+             !
+             !-- USM surfaces
+             IF ( surf_usm_h%ns > 0 )  THEN
 
 
@@ -1604,37 +1504,39 @@
                    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
+                      !
+                      !-- 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)*        & 
-                         !                                              kg/m^3
-                                                                       rho_air(k)   
+                         !
+                         !--          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
 
  
                       ELSE
-
-                         !> VOCs
-                         IF ( len_index_voc .GT. 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) *   &
-                                                                          !    m^3/mole          ppm
-                                                                          conv_to_ratio(k,j,i)*ratio2ppm*    &
-                         !                                                kg/m^3
-                                                                          rho_air(k)   
-
-                         !> OTHER SPECIES
+                         
+                         !
+                         !-- 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
                          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) *      &
-                                                                             !  mole/Kg    
-                                                                         (1/emt_att%xm(ispec))*                 &
-                            !                                                m^3/mole           ppm
-                                                                         conv_to_ratio(k,j,i)*ratio2ppm*        &
-                            !                                            kg/m^3
-                                                                         rho_air(k)   
+                            !
+                            !--          (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   
 
 
@@ -1646,16 +1548,16 @@
  
                 ENDDO
-             END IF
+
+             ENDIF
+
           ENDDO
 
-       ENDIF 
-
-
-    !> At the end of each call to chem_emissions setup, the time_factor is deallocated (ALLOCATED ONLY in the DEFAULT mode) 
-
+       ENDIF
+
+       !
+       !-- Deallocate time_factor in case of DEFAULT mode) 
        IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor )
 
-   ENDIF !> emis_output_required
-
+   ENDIF
 
  END SUBROUTINE chem_emissions_setup
Index: palm/trunk/SOURCE/chem_modules.f90
===================================================================
--- palm/trunk/SOURCE/chem_modules.f90	(revision 3610)
+++ palm/trunk/SOURCE/chem_modules.f90	(revision 3611)
@@ -27,7 +27,10 @@
 ! -----------------
 ! $Id$
+! Minor formatting
+!
+! 3458 2018-10-30 14:51:23Z kanani
 ! from chemistry branch r3443:
 ! ??
-! 
+!
 ! 3298 2018-10-02 12:21:11Z kanani
 ! - Minor formatting
@@ -77,68 +80,69 @@
     PUBLIC spc_names
 
-    INTEGER(iwp), DIMENSION(99)  :: cs_pr_index            = 0
-    INTEGER(iwp)  :: ibc_cs_b                                                      !< integer flag for bc_cs_b
-    INTEGER(iwp)  :: ibc_cs_t                                                      !< integer flag for bc_cs_t
-    INTEGER(iwp)  :: cs_pr_count                           = 0 
-    INTEGER(iwp)  :: max_pr_cs                             = 0
-    INTEGER(iwp)  :: cs_vertical_gradient_level_ind(99,10) = -9999                 !< grid index values of cs_vertical_gradient_level_ind(s)
-
-    LOGICAL       :: constant_top_csflux(99)               = .TRUE.                !< chem spcs at the top  orig .TRUE.
-    LOGICAL       :: constant_csflux(99)                   = .TRUE.                !< chem spcs at namelist parameter   orig TRUE
-    LOGICAL       :: call_chem_at_all_substeps             = .FALSE.               !< namelist parameter 
-    LOGICAL       :: chem_debug0                           = .FALSE.               !< namelist parameter flag for minimum print output
-    LOGICAL       :: chem_debug1                           = .FALSE.               !< namelist parameter flag for print output
-    LOGICAL       :: chem_debug2                           = .FALSE.               !< namelist parameter flag for further print output
-    LOGICAL       :: chem_gasphase_on                      = .TRUE.                !< namelist parameter 
-    LOGICAL       :: emission_output_required              = .TRUE.                !< Logical Variable for requiring Emission Outputs
-    LOGICAL       :: do_emis                               = .FALSE.               !< Flag for turning on chemistry emissions 
-    LOGICAL       :: cs_pr_namelist_found                  = .FALSE.               !< Namelist parameter: Names of t
-    LOGICAL       :: do_depo                               = .FALSE.               !< namelist parameter for activation of deposition calculation
-
-
-
+    INTEGER(iwp), DIMENSION(99) :: cs_pr_index            = 0
+    INTEGER(iwp) ::  ibc_cs_b                                                      !< integer flag for bc_cs_b
+    INTEGER(iwp) ::  ibc_cs_t                                                      !< integer flag for bc_cs_t
+    INTEGER(iwp) ::  cs_pr_count                           = 0 
+    INTEGER(iwp) ::  max_pr_cs                             = 0
+    INTEGER(iwp) ::  cs_vertical_gradient_level_ind(99,10) = -9999                 !< grid index values of cs_vertical_gradient_level_ind(s)
+
+    LOGICAL      ::  constant_top_csflux(99)               = .TRUE.                !< chem spcs at the top  orig .TRUE.
+    LOGICAL      ::  constant_csflux(99)                   = .TRUE.                !< chem spcs at namelist parameter   orig TRUE
+    LOGICAL      ::  call_chem_at_all_substeps             = .FALSE.               !< namelist parameter 
+    LOGICAL      ::  chem_debug0                           = .FALSE.               !< namelist parameter flag for minimum print output
+    LOGICAL      ::  chem_debug1                           = .FALSE.               !< namelist parameter flag for print output
+    LOGICAL      ::  chem_debug2                           = .FALSE.               !< namelist parameter flag for further print output
+    LOGICAL      ::  chem_gasphase_on                      = .TRUE.                !< namelist parameter 
+    LOGICAL      ::  emission_output_required              = .TRUE.                !< Logical Variable for requiring Emission Outputs
+    LOGICAL      ::  do_emis                               = .FALSE.               !< Flag for turning on chemistry emissions 
+    LOGICAL      ::  cs_pr_namelist_found                  = .FALSE.               !< Namelist parameter: Names of t
+    LOGICAL      ::  do_depo                               = .FALSE.               !< namelist parameter for activation of deposition calculation
+
+
+!
 !-- Namelist parameters for creating initial chemistry profiles
-    REAL(wp) :: wall_csflux (99,0:5)               = 0.0_wp                        !< namelist parameter
-    REAL(wp) :: cs_vertical_gradient (99,10)       = 0.0_wp                        !< namelist parameter
-    REAL(wp) :: cs_vertical_gradient_level (99,10) = -999999.9_wp                  !< namelist parameter
-    REAL(wp) :: top_csflux ( 99 )                  = 0.0_wp                        !< namelist parameter 
-    REAL(wp) :: cs_surface_initial_change(99)      = 0.0_wp                        !< namelist parameter
-    REAL(wp) :: surface_csflux(99 )                = 0.0_wp                        !< namelist parameter: fluxes where 'surface_csflux_name' is in the namelist
-
-    REAL(wp), DIMENSION(:),  ALLOCATABLE              :: bc_cs_t_val
-    REAL(wp), DIMENSION(:),  ALLOCATABLE              ::  css                      !< scaling parameter for chem spcs
-    REAL(wp), DIMENSION(99)                           :: cs_surface = 0.0_wp       !< Namelist parameter: Surface conc of chem spcs'
-    REAL(wp), DIMENSION(99,100)                       :: cs_heights = 9999999.9_wp !< Namelist parameter: Height lvls(m) for cs_profiles
-    REAL(wp), DIMENSION(99,100)                       :: cs_profile = 9999999.9_wp !< Namelist parameter: Chem conc for each spcs defined
+    REAL(wp) ::  wall_csflux (99,0:5)               = 0.0_wp                        !< namelist parameter
+    REAL(wp) ::  cs_vertical_gradient (99,10)       = 0.0_wp                        !< namelist parameter
+    REAL(wp) ::  cs_vertical_gradient_level (99,10) = -999999.9_wp                  !< namelist parameter
+    REAL(wp) ::  top_csflux ( 99 )                  = 0.0_wp                        !< namelist parameter 
+    REAL(wp) ::  cs_surface_initial_change(99)      = 0.0_wp                        !< namelist parameter
+    REAL(wp) ::  surface_csflux(99 )                = 0.0_wp                        !< namelist parameter: fluxes where 'surface_csflux_name' is in the namelist
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE               ::  bc_cs_t_val
+    REAL(wp), DIMENSION(:), ALLOCATABLE               ::  css                       !< scaling parameter for chem spcs
+    REAL(wp), DIMENSION(99)                           ::  cs_surface = 0.0_wp       !< Namelist parameter: Surface conc of chem spcs'
+    REAL(wp), DIMENSION(99,100)                       ::  cs_heights = 9999999.9_wp !< Namelist parameter: Height lvls(m) for cs_profiles
+    REAL(wp), DIMENSION(99,100)                       ::  cs_profile = 9999999.9_wp !< Namelist parameter: Chem conc for each spcs defined
 
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET   :: cs                        !< chem spcs
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET   :: cs_p                      !< prognostic value of chem spc
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET   :: tcs_m                     !< weighted tendency of cs for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET   ::  cs                        !< chem spcs
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET   ::  cs_p                      !< prognostic value of chem spc
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET   ::  tcs_m                     !< weighted tendency of cs for previous sub-timestep (Runge-Kutta)
 
 #else                                                               
-! use pointers cs, cs_p and tcs_m to point arrays cs_1, cs_2, and cs_3
-    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: cs_1                      !< pointer for swapping of timelevels for respective quantity
-    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: cs_2                      !< pointer for swapping of timelevels for respective quantity
-    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: cs_3                      !< pointer for swapping of timelevels for respective quantity
-    REAL(wp), DIMENSION(:,:,:), POINTER               :: cs                        !< pointer: sgs chem spcs)
-    REAL(wp), DIMENSION(:,:,:), POINTER               :: cs_p                      !< pointer: prognostic value of sgs chem spcs 
-    REAL(wp), DIMENSION(:,:,:), POINTER               :: tcs_m                     !< pointer:
+!
+!-- Use pointers cs, cs_p and tcs_m to point arrays cs_1, cs_2, and cs_3
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  cs_1                      !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  cs_2                      !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  cs_3                      !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), POINTER               ::  cs                        !< pointer: sgs chem spcs)
+    REAL(wp), DIMENSION(:,:,:), POINTER               ::  cs_p                      !< pointer: prognostic value of sgs chem spcs 
+    REAL(wp), DIMENSION(:,:,:), POINTER               ::  tcs_m                     !< pointer:
 
 #endif                                                                            
  
-    CHARACTER (LEN=20)                 :: bc_cs_b             = 'dirichlet'        !< namelist parameter
-    CHARACTER (LEN=20)                 :: bc_cs_t             = 'initial_gradient' !< namelist parameter
-    CHARACTER (LEN=11), DIMENSION(99)  :: cs_name             = 'novalue'          !< Namelist parameter: chem spcs names 
-    CHARACTER (LEN=11), DIMENSION(99)  :: cs_profile_name     = 'novalue'          !< Namelist parameter: Names of the chem for profiles
-    CHARACTER (LEN=11), DIMENSION(99)  :: surface_csflux_name = 'novalue'          !< Namelist parameter: chem species surface fluxes names
-                                                                                   !< active chem spcs, default is 'novalue')  ????
-    CHARACTER (LEN=80)                 :: mode_emis           ='PARAMETERIZED'     !< Mode of chemistry emissions: DEFAULT .OR. EXPERT .OR. 
-                                                                                   !  PARAMETERIZED
-    CHARACTER (LEN=80)                 :: time_fac_type       ='MDH'               !< Type of time treatment in the emis DEFAULT mode: HOUR .OR. MDH 
-    CHARACTER (LEN=80)                 :: daytype_mdh         ='workday'           !< Type of day in the MDH case: workday, weekend, holiday
-    CHARACTER (LEN=11), DIMENSION(99)  :: data_output_pr_cs   = 'novalue'          !< Namelist parameter: Names of the che    m for profile output
-                                                                                   !< by cs_name for each height lvls defined by cs_heights
+    CHARACTER (LEN=20)                ::  bc_cs_b             = 'dirichlet'         !< namelist parameter
+    CHARACTER (LEN=20)                ::  bc_cs_t             = 'initial_gradient'  !< namelist parameter
+    CHARACTER (LEN=11), DIMENSION(99) ::  cs_name             = 'novalue'           !< Namelist parameter: chem spcs names 
+    CHARACTER (LEN=11), DIMENSION(99) ::  cs_profile_name     = 'novalue'           !< Namelist parameter: Names of the chem for profiles
+    CHARACTER (LEN=11), DIMENSION(99) ::  surface_csflux_name = 'novalue'           !< Namelist parameter: chem species surface fluxes names
+                                                                                    !< active chem spcs, default is 'novalue')  ????
+    CHARACTER (LEN=80)                ::  mode_emis           ='PARAMETERIZED'      !< Mode of chemistry emissions: DEFAULT .OR. EXPERT .OR. 
+                                                                                    !< PARAMETERIZED
+    CHARACTER (LEN=80)                ::  time_fac_type       ='MDH'                !< Type of time treatment in the emis DEFAULT mode: HOUR .OR. MDH 
+    CHARACTER (LEN=80)                ::  daytype_mdh         ='workday'            !< Type of day in the MDH case: workday, weekend, holiday
+    CHARACTER (LEN=11), DIMENSION(99) ::  data_output_pr_cs   = 'novalue'           !< Namelist parameter: Names of the che    m for profile output
+                                                                                    !< by cs_name for each height lvls defined by cs_heights
 !
 !-- Namelist parameters for chem_emissions
@@ -150,54 +154,52 @@
     REAL(wp) ::  emiss_factor_main ( 99 ) = -9999.0_wp
     REAL(wp) ::  emiss_factor_side ( 99 ) = -9999.0_wp
-    
+!    
 !-- Other Emissions Variables
     INTEGER(iwp) ::  nspec_out                                                     !< Output of routine chem_emis_matching with 
                                                                                    !< number of matched species
-    REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:)         :: emis_distribution          !> Emissions Final Values (main module output)
-
-    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)            :: match_spec_input           !< Index of Input chem species for matching routine 
-    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)            :: match_spec_model           !< Index of Model chem species for matching routine
-    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)            :: match_spec_voc_input       !< index of VOC input components matching the model's VOCs
-    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)            :: match_spec_voc_model       !< index of VOC model species matching the input VOCs comp.
-    INTEGER(iwp),DIMENSION(:)                        :: match_spec_pm(1:3)         !< results of matching the input and model's PMs
-    INTEGER(iwp),DIMENSION(:)                        :: match_spec_nox(1:2)        !< results of matching the input and model's NOx
-    INTEGER(iwp),DIMENSION(:)                        :: match_spec_sox(1:2)        !< results of matching the input and model's SOx! 
-                                                                                   ! TBD: evaluate whether to make them allocatable 
-                                                                                   ! and allocate their dimension in the matching 
-                                                                                   ! routine according to the number of components 
-                                                                                   ! matching between the model and the input files
-
+    REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:) ::  emis_distribution                 !> Emissions Final Values (main module output)
+
+    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)    ::  match_spec_input                  !< Index of Input chem species for matching routine 
+    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)    ::  match_spec_model                  !< Index of Model chem species for matching routine
+    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)    ::  match_spec_voc_input              !< index of VOC input components matching the model's VOCs
+    INTEGER(iwp),ALLOCATABLE,DIMENSION(:)    ::  match_spec_voc_model              !< index of VOC model species matching the input VOCs comp.
+    INTEGER(iwp),DIMENSION(:)                ::  match_spec_pm(1:3)                !< results of matching the input and model's PMs
+    INTEGER(iwp),DIMENSION(:)                ::  match_spec_nox(1:2)               !< results of matching the input and model's NOx
+    INTEGER(iwp),DIMENSION(:)                ::  match_spec_sox(1:2)               !< results of matching the input and model's SOx! 
+                                                                                  
+
+!
+!-- Selected atomic/molecular weights:
+    REAL, PARAMETER        ::  xm_H     =    1.00790e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_N     =   14.00670e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_C     =   12.01115e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_S     =   32.06400e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_O     =   15.99940e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_F     =   18.99840e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Na    =   22.98977e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Cl    =   35.45300e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Rn222 =  222.00000e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Pb210 =  210.00000e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Ca    =   40.07800e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_K     =   39.09800e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Mg    =   24.30500e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Pb    =  207.20000e-3           !< kg/mol
+    REAL, PARAMETER        ::  xm_Cd    =  112.41000e-3           !< kg/mol
     
-    !-- Selected atomic/molecular weights:
+    REAL, PARAMETER        ::  xm_h2o   = xm_H * 2 + xm_O         !< kg/mol
+    REAL, PARAMETER        ::  xm_o3    = xm_O * 3                !< kg/mol
+    REAL, PARAMETER        ::  xm_N2O5  = xm_N * 2 + xm_O * 5     !< kg/mol
+    REAL, PARAMETER        ::  xm_HNO3  = xm_H + xm_N + xm_O * 3  !< kg/mol
+    REAL, PARAMETER        ::  xm_NH4   = xm_N + xm_H * 4         !< kg/mol
+    REAL, PARAMETER        ::  xm_SO4   = xm_S + xm_O * 4         !< kg/mol
+    REAL, PARAMETER        ::  xm_NO3   = xm_N + xm_O * 3         !< kg/mol
+    REAL, PARAMETER        ::  xm_CO2   = xm_C + xm_O * 2         !< kg/mol
     
-    REAL, PARAMETER        ::  xm_H     =    1.00790e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_N     =   14.00670e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_C     =   12.01115e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_S     =   32.06400e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_O     =   15.99940e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_F     =   18.99840e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Na    =   22.98977e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Cl    =   35.45300e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Rn222 =  222.00000e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Pb210 =  210.00000e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Ca    =   40.07800e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_K     =   39.09800e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Mg    =   24.30500e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Pb    =  207.20000e-3           ! kg/mol
-    REAL, PARAMETER        ::  xm_Cd    =  112.41000e-3           ! kg/mol
-    
-    REAL, PARAMETER        ::  xm_h2o   = xm_H * 2 + xm_O         ! kg/mol
-    REAL, PARAMETER        ::  xm_o3    = xm_O * 3                ! kg/mol
-    REAL, PARAMETER        ::  xm_N2O5  = xm_N * 2 + xm_O * 5     ! kg/mol
-    REAL, PARAMETER        ::  xm_HNO3  = xm_H + xm_N + xm_O * 3  ! kg/mol
-    REAL, PARAMETER        ::  xm_NH4   = xm_N + xm_H * 4         ! kg/mol
-    REAL, PARAMETER        ::  xm_SO4   = xm_S + xm_O * 4         ! kg/mol
-    REAL, PARAMETER        ::  xm_NO3   = xm_N + xm_O * 3         ! kg/mol
-    REAL, PARAMETER        ::  xm_CO2   = xm_C + xm_O * 2         ! kg/mol
-    
-    ! mass of air
-    REAL, PARAMETER        ::  xm_air   =  28.964e-3              ! kg/mol
+!
+!-- mass of air
+    REAL, PARAMETER        ::  xm_air   =  28.964e-3              !< kg/mol
         
-    ! dummy weight, used for complex molecules:
+!
+!-- dummy weight, used for complex molecules:
     REAL, PARAMETER        ::  xm_dummy =  1000.0e-3              ! kg/mol
 
Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3610)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3611)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Minor formatting             
+!
+! 3600 2018-12-04 13:49:07Z banzhafs
 ! Code update to comply PALM coding rules            
 ! Bug fix in par_dir_diff subroutine                 
@@ -825,5 +828,5 @@
     INTEGER(iwp) ::  ss        !< 
     REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
-    REAL(wp), DIMENSION(nzb:nzt+1, nysg:nyng, nxlg:nxrg) ::  cs_3d
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  cs_3d
     REAL(wp) ::  flag !< flag to mask topography grid points
 
@@ -1207,5 +1210,5 @@
     LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
     REAL(wp)     ::  fill_value 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf  
 
     !
@@ -1294,5 +1297,5 @@
 
     REAL(wp)             ::  fill_value   !<
-    REAL(sp), DIMENSION(nxl:nxr, nys:nyn, nzb_do:nzt_do) ::  local_pf 
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
 
 
@@ -2359,9 +2362,9 @@
 
     INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
-    REAL(wp), DIMENSION(nzb+1:nzt, 0:threads_per_task-1)          :: flux_s_cs   !<
-    REAL(wp), DIMENSION(nzb+1:nzt, 0:threads_per_task-1)          :: diss_s_cs   !<
-    REAL(wp), DIMENSION(nzb+1:nzt, nys:nyn, 0:threads_per_task-1) :: flux_l_cs   !<
-    REAL(wp), DIMENSION(nzb+1:nzt, nys:nyn, 0:threads_per_task-1) :: diss_l_cs   !<
-    REAL(wp), DIMENSION(0:nz+1)                                 :: pr_init_cs  !<
+    REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)          :: flux_s_cs   !<
+    REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)          :: diss_s_cs   !<
+    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1)  :: flux_l_cs   !<
+    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1)  :: diss_l_cs   !<
+    REAL(wp), DIMENSION(0:nz+1)                                  :: pr_init_cs  !<
 
     !
@@ -2473,5 +2476,5 @@
     LOGICAL, INTENT(OUT) :: found 
 
-    REAL(wp), DIMENSION(nzb:nzt+1, nys_on_file-nbgp:nyn_on_file+nbgp, nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !< 3D array to temp store data
+    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !< 3D array to temp store data
 
 
