!> @file chem_emissions_mod.f90
!--------------------------------------------------------------------------------!
! This file is part of PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 2018-2018 Leibniz Universitaet Hannover
! Copyright 2018-2018 Freie Universitaet Berlin
! Copyright 2018-2018 Karlsruhe Institute of Technology
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: chem_emissions_mod.f90 3772 2019-02-28 15:51:57Z forkel $
! - In case of parametrized emissions, assure that emissions are only on natural
!   surfaces (i.e. streets) and not on urban surfaces. 
! - some unnecessary if clauses removed
! 
! 3685 2019-01-21 01:02:11Z knoop
! Some interface calls moved to module_interface + cleanup
! 
! 3611 2018-12-07 14:14:11Z banzhafs
! 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:
! Break lines at 132 characters
!
! 3483 2018-11-02 14:19:26Z raasch
! bugfix: wrong locations of netCDF directives fixed
! 
! 3467 2018-10-30 19:05:21Z suehring
! Enabled PARAMETRIZED mode for default surfaces when LSM is not applied but 
! salsa is used
! 
! 3458 2018-10-30 14:51:23Z kanani
! from chemistry branch r3443, banzhafs, Russo:
! Additional correction for index of input file of pre-processed mode
! Removed atomic and molecular weights as now available in chem_modules and
! added par_emis_time_factor (formerly in netcdf_data_input_mod)
! Added correction for index of input file of pre-processed mode
! Added a condition for default mode necessary for information read from ncdf file
! in pre-processed and default mode
! Correction of multiplying factor necessary for scaling emission values in time
! Introduced correction for scaling units in the case of DEFAULT emission mode
! 
! 3373 2018-10-18 15:25:56Z kanani
! Fix wrong location of __netcdf directive
! 
! 3337 2018-10-12 15:17:09Z kanani
! (from branch resler)
! Formatting
! 
! 3312 2018-10-06 14:15:46Z knoop
! Initial revision
!
! 3286 2018-09-28 07:41:39Z forkel
!
! Authors:
! --------
! @author Emmanuele Russo (FU-Berlin)
! @author Sabine Banzhaf  (FU-Berlin)
! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
! 
! Description:
! ------------
!>  MODULE for reading-in Chemistry Emissions data
!>
!> @todo Check_parameters routine should be developed: for now it includes just one condition
!> @todo Use of Restart files not contempled at the moment
!> @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>
!------------------------------------------------------------------------------!

 MODULE chem_emissions_mod

    USE arrays_3d,                                                             &
        ONLY:  rho_air

    USE control_parameters,                                                    &
        ONLY:  end_time, message_string, initializing_actions,                 &
               intermediate_timestep_count, dt_3d
 
    USE indices

    USE kinds

#if defined( __netcdf )
    USE netcdf
#endif

    USE netcdf_data_input_mod,                                                  &
        ONLY: chem_emis_att_type, chem_emis_val_type

    USE date_and_time_mod,                                                      &
        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
 
    USE chem_modules

    USE statistics,                                                             &   
        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
                           
    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


!    
!-- Checks Input parameters
    INTERFACE chem_emissions_check_parameters
       MODULE PROCEDURE chem_emissions_check_parameters
    END INTERFACE chem_emissions_check_parameters
!
!-- 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
       MODULE PROCEDURE chem_emissions_setup
    END INTERFACE chem_emissions_setup


    
    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
!
!-- Public Variables
    PUBLIC con_factor, len_index, len_index_pm, len_index_voc

 CONTAINS

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Routine for checking input parameters
!------------------------------------------------------------------------------!
 SUBROUTINE chem_emissions_check_parameters


    IMPLICIT NONE

    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
    

 END SUBROUTINE chem_emissions_check_parameters

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Matching the chemical species indices. The routine checks what are the indices of the emission input species 
!> and the corresponding ones of the model species. The routine gives as output a vector containing the number 
!> of common species: it is important to note that while the model species are distinct, their values could be 
!> given to a single species in input: for example, in the case of NO2 and NO, values may be passed in input as
!> NOx values.
!------------------------------------------------------------------------------!
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


    CALL location_message( 'Matching input emissions and model chemistry species', .FALSE. )

    !
    !-- Number of input emission species.
    nspec_emis_inp=emt_att%nspec 

    !
    !-- Check the emission mode: DEFAULT, PRE-PROCESSED or PARAMETERIZED 
    SELECT CASE( TRIM( mode_emis ) ) 

       !
       !-- PRE-PROCESSED mode
       CASE ( "PRE-PROCESSED" )

          len_index = 0
          len_index_voc = 0

          IF ( nvar > 0 .AND. (nspec_emis_inp > 0) )  THEN
             !
             !-- Cycle over model species
             DO ind_mod = 1,  nvar 
                !
                !-- Cycle over input species  
                DO ind_inp = 1, nspec_emis_inp

                   !
                   !-- Check for VOC Species  
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN        
                      DO ind_voc = 1, emt_att%nvoc
             
                         IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
                            len_index = len_index + 1
                            len_index_voc = len_index_voc + 1
                         ENDIF
                      END DO
                   ENDIF
                   !
                   !-- Other Species
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                      len_index = len_index + 1
                   ENDIF
                ENDDO
             ENDDO

             !
             !-- Allocate array for storing the indices of the matched species
             IF ( len_index > 0 )  THEN
 
                ALLOCATE( match_spec_input(len_index) ) 
 
                ALLOCATE( match_spec_model(len_index) )

                IF ( len_index_voc > 0 )  THEN
                   !
                   !-- contains indices of the VOC model species
                   ALLOCATE( match_spec_voc_model(len_index_voc) )
                   !
                   !-- contains the indices of different values of VOC composition of input variable VOC_composition
                   ALLOCATE( match_spec_voc_input(len_index_voc) )

                ENDIF

                !
                !-- pass the species indices to declared arrays
                len_index = 0

                !
                !-- Cycle over model species
                DO ind_mod = 1, nvar 
                   !
                   !-- Cycle over Input species  
                   DO ind_inp = 1, nspec_emis_inp
                      !
                      !-- VOCs
                      IF ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND.    &
                           ALLOCATED(match_spec_voc_input) )  THEN
                         
                         DO ind_voc= 1, emt_att%nvoc
                            IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
                               len_index = len_index + 1
                               len_index_voc = len_index_voc + 1
                        
                               match_spec_input(len_index) = ind_inp
                               match_spec_model(len_index) = ind_mod

                               match_spec_voc_input(len_index_voc) = ind_voc
                               match_spec_voc_model(len_index_voc) = ind_mod                         
                            ENDIF
                         END DO
                      ENDIF

                      !
                      !-- Other Species
                      IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                         len_index = len_index + 1
                         match_spec_input(len_index) = ind_inp
                         match_spec_model(len_index) = ind_mod
                      ENDIF
                   END DO
                END DO

             ELSE
                !
                !-- in case there are no species matching
                message_string = 'Non of given emission species'            //         &
                                 ' matches'                                //          &
                                 ' model chemical species:'                //          &
                                 ' Emission routine is not called'  
                CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
             ENDIF
 
          ELSE

             !
             !-- either spc_names is zero or nspec_emis_inp is not allocated
             message_string = 'Array of Emission species not allocated:'             //          &
                              ' Either no emission species are provided as input or'  //         &
                              ' no chemical species are used by PALM:'                //         &
                              ' Emission routine is not called'                  
             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) 

          ENDIF

       !
       !-- DEFAULT mode
       CASE ("DEFAULT")

          len_index = 0          !<  index for TOTAL number of species   
          len_index_voc = 0      !<  index for TOTAL number of VOCs
          len_index_pm = 3       !<  index for TOTAL number of PMs: PM1, PM2.5, PM10.
 
          IF ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN

             !
             !-- Cycle over model species
             DO ind_mod = 1, nvar
                !
                !-- Cycle over input species
                DO ind_inp = 1, nspec_emis_inp

                   !
                   !-- Check for VOC Species  
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
                      DO ind_voc= 1, emt_att%nvoc
                        
                         IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
                            len_index = len_index + 1
                            len_index_voc = len_index_voc + 1
                         ENDIF
                         
                      END DO
                   ENDIF

                   !
                   !-- PMs: There is one input species name for all PM 
                   !-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
                      !
                      !-- PM1
                      IF ( TRIM( spc_names(ind_mod) ) == "PM1" )  THEN
                         len_index = len_index + 1
                      !
                      !-- PM2.5
                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
                         len_index = len_index + 1
                      !
                      !-- PM10
                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
                         len_index = len_index + 1
                      ENDIF
                   ENDIF

                   !
                   !-- NOx: NO2 and NO    
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" )  THEN 
                      !
                      !-- NO
                      IF ( TRIM( spc_names(ind_mod) ) == "NO" )  THEN
                         len_index = len_index + 1
                      !
                      !-- NO2
                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
                         len_index = len_index + 1
                      ENDIF
                   ENDIF

                   !
                   !-- SOx: SO2 and SO4
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" )  THEN
                      !
                      !-- SO2
                      IF ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
                         len_index = len_index + 1
                      !
                      !-- SO4
                      ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
                         len_index = len_index + 1
                      ENDIF
                   ENDIF

                   !
                   !-- Other Species
                   IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                      len_index = len_index + 1
                   ENDIF
                END DO
             END DO


             !
             !-- Allocate arrays
             IF ( len_index > 0 )  THEN

                ALLOCATE( match_spec_input(len_index) )
                ALLOCATE( match_spec_model(len_index) )

                IF ( len_index_voc > 0 )  THEN
                   !
                   !-- Contains indices of the VOC model species
                   ALLOCATE( match_spec_voc_model(len_index_voc) )
                   !
                   !-- Contains the indices of different values of VOC composition
                   !-- of input variable VOC_composition
                   ALLOCATE( match_spec_voc_input(len_index_voc) )                                                  
                ENDIF

                !
                !-- Pass the species indices to declared arrays
                len_index = 0
                len_index_voc = 0
                
                DO ind_mod = 1, nvar
                   DO ind_inp = 1, nspec_emis_inp 

                      !
                      !-- VOCs
                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND.   &
                         ALLOCATED(match_spec_voc_input) )  THEN      
                         DO ind_voc= 1, emt_att%nvoc
                            IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
                               len_index = len_index + 1
                               len_index_voc = len_index_voc + 1
                        
                               match_spec_input(len_index) = ind_inp
                               match_spec_model(len_index) = ind_mod

                               match_spec_voc_input(len_index_voc) = ind_voc
                               match_spec_voc_model(len_index_voc) = ind_mod                         
                            ENDIF
                         END DO
                      ENDIF

                      !
                      !-- PMs 
                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
                         !
                         !-- PM1
                         IF ( TRIM( spc_names(ind_mod) ) == "PM1" )  THEN
                            len_index = len_index + 1

                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
                         !
                         !-- PM2.5
                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" )  THEN
                            len_index = len_index + 1

                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
                         !
                         !-- PM10
                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" )  THEN
                            len_index = len_index + 1
    
                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
  
                         ENDIF
                      ENDIF

                      !
                      !-- NOx 
                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOx" )  THEN
                         !
                         !-- NO
                         IF ( TRIM( spc_names(ind_mod) ) == "NO" )  THEN
                            len_index = len_index + 1

                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
                         !
                         !-- NO2
                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" )  THEN
                            len_index = len_index + 1

                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
  
                         ENDIF
                      ENDIF

                      !
                      !-- SOx 
                      IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOx" ) THEN
                         !
                         !-- SO2
                         IF ( TRIM( spc_names(ind_mod) ) == "SO2" )  THEN
                            len_index = len_index + 1

                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
  
                         !
                         !-- SO4
                         ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" )  THEN
                            len_index = len_index + 1

                            match_spec_input(len_index) = ind_inp
                            match_spec_model(len_index) = ind_mod
  
                         ENDIF
                      ENDIF

                      !
                      !-- Other Species
                      IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                         len_index = len_index + 1
                            
                         match_spec_input(len_index) = ind_inp
                         match_spec_model(len_index) = ind_mod
                      ENDIF
                   END DO
                END DO

             ELSE

                message_string = 'Non of given Emission Species'            //         &
                                 ' matches'                                //          &
                                 ' model chemical species'                 //          &
                                 ' Emission routine is not called'         
                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 

             ENDIF

          ELSE

             message_string = 'Array of Emission species not allocated: '            //          &
                              ' Either no emission species are provided as input or'  //         &
                              ' no chemical species are used by PALM:'                //         &
                              ' Emission routine is not called'                                   
             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) 
 
          ENDIF
 
       !
       !-- PARAMETERIZED mode
       CASE ("PARAMETERIZED")

          len_index = 0

          IF ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN

             !
             !-- Cycle over model species
             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
                ENDDO
             ENDDO

             IF ( len_index > 0 )  THEN

                !
                !-- Allocation of Arrays of the matched species
                ALLOCATE( match_spec_input(len_index) ) 
 
                ALLOCATE( match_spec_model(len_index) )

                !
                !-- Pass the species indices to declared arrays    
                len_index = 0

                DO ind_mod = 1, nvar  
                   ind_inp = 1
                   DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )     
                      IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                         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
                   END DO
                END DO

                !
                !-- 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'    //          &
                                       ' emiss_factor_main AND emiss_factor_side'           //          &
                                       ' not provided for each of the emissions species'    //          &
                                       ' or not provided at all: PLEASE set a finite value' //          &
                                       ' for these parameters in the chemistry namelist'          
                      CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 ) 
                   ENDIF
                END DO


             ELSE
                
                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

          ELSE
     
             message_string = 'Array of Emission species not allocated: '            //           &
                              ' Either no emission species are provided as input or'  //          &
                              ' no chemical species are used by PALM.'                //          &
                              ' Emission routine is not called'                                   
             CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) 
 
          ENDIF             


       !
       !-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
       CASE DEFAULT

          message_string = 'Emission Module switched ON, but'            //                         &
                           ' either no emission mode specified or incorrectly given :'  //          &
                           ' please, pass the correct value to the namelist parameter "mode_emis"'                  
          CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )              

       END SELECT

 END SUBROUTINE chem_emissions_match

 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization:
!> Netcdf reading, arrays allocation and first calculation of cssws
!> fluxes at timestep 0 
!------------------------------------------------------------------------------!
 SUBROUTINE chem_emissions_init

    USE surface_mod,                                                           &
        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h

    USE netcdf_data_input_mod,                                                 &
        ONLY:  chem_emis, chem_emis_att
    
    IMPLICIT NONE
 
    CHARACTER (LEN=80)          ::  units                                           !< units of inputs

    INTEGER(iwp)                :: ispec                                            !< running index

!    
!-- 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( chem_emis_att, nspec_out )
 
  IF ( nspec_out == 0 )  THEN
  
     emission_output_required = .FALSE.

  ELSE

     emission_output_required = .TRUE.


     !
     !-- Set molecule masses' 
     ALLOCATE( chem_emis_att%xm(nspec_out) )

     DO ispec = 1, nspec_out
        SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
           CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2        !< kg/mole
           CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4        !< kg/mole
           CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O             !< kg/mole
           CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2        !< kg/mole   
           CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3        !< kg/mole
           CASE ( 'CO'  ); chem_emis_att%xm(ispec) = xm_C + xm_O            !< kg/mole
           CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2        !< kg/mole
           CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4        !< kg/mole     
           CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3  !< kg/mole  
           CASE DEFAULT
              chem_emis_att%xm(ispec) = 1.0_wp
        END SELECT
     ENDDO

    
     !
     !-- assign emission values
     SELECT CASE ( TRIM( mode_emis ) )   


        !
        !-- PRE-PROCESSED case
        CASE ( "PRE-PROCESSED" )

           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,nspec_out) )  
 
           !
           !-- Get emissions at the first time step
           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )

        !
        !-- Default case
        CASE ( "DEFAULT" )

           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) ) 

           !
           !-- Get emissions at the first time step
           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )

        !
        !-- PARAMETERIZED case
        CASE ( "PARAMETERIZED" )

           IF ( .NOT. ALLOCATED( emis_distribution) ) ALLOCATE( emis_distribution(1,0:ny,0:nx,nspec_out) )

           !
           !-- Get emissions at the first time step
           CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out)

     END SELECT

  ENDIF    

 END SUBROUTINE chem_emissions_init



!------------------------------------------------------------------------------!
! Description:
! ------------
!> Routine for Update of Emission values at each timestep 
!-------------------------------------------------------------------------------!

 SUBROUTINE chem_emissions_setup( emt_att, emt, nspec_out )

   USE arrays_3d,                                                    &
       ONLY:  dzw
   USE grid_variables,                                               &
       ONLY: dx, dy
   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,                                                    &        
       ONLY: hyp, pt 

   
 IMPLICIT NONE
  

    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
    
    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  
       IF ( call_chem_at_all_substeps )  THEN

          dt_emis = dt_3d * weight_pres(intermediate_timestep_count)

       ELSE

          dt_emis = dt_3d

       ENDIF


       !
       !-- Conversion of units to the ones employed in PALM 
       !-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed

        IF ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN

          SELECT CASE ( TRIM( emt_att%units ) )
             !
             !-- kilograms
             CASE ( 'kg/m2/s', 'KG/M2/S' ) 
                
                con_factor=1

             CASE ('kg/m2/hour', 'KG/M2/HOUR' )

                con_factor=hour_to_s

             CASE ( 'kg/m2/day', 'KG/M2/DAY' ) 
                
                con_factor=day_to_s

             CASE ( 'kg/m2/year', 'KG/M2/YEAR' )
              
                con_factor=year_to_s

             !
             !-- Tons
             CASE ( 'ton/m2/s', 'TON/M2/S' ) 
               
                con_factor=tons_to_kg

             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ) 
               
                con_factor=tons_to_kg*hour_to_s

             CASE ( 'ton/m2/year', 'TON/M2/YEAR' ) 
               
                con_factor=tons_to_kg*year_to_s

             !
             !-- Grams
             CASE ( 'g/m2/s', 'G/M2/S' )
          
                con_factor=g_to_kg

             CASE ( 'g/m2/hour', 'G/M2/HOUR' ) 

                con_factor=g_to_kg*hour_to_s

             CASE ( 'g/m2/year', 'G/M2/YEAR' )
 
                con_factor=g_to_kg*year_to_s

             !
             !-- Micrograms
             CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' )
 
                con_factor=miug_to_kg

             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ) 
 
                con_factor=miug_to_kg*hour_to_s

             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' )
  
                con_factor=miug_to_kg*year_to_s

             CASE DEFAULT  
                message_string = 'The Units of the provided emission input' // &
                                 ' are not the ones required by PALM-4U: please check '      // &
                                 ' emission module documentation.'                                 
                CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )

          END SELECT

         
       ENDIF

       !
       !-- 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 multiply the input by (m**2/N)*dz
             !>  (m**2/N)*dz == V/N 
             !>  V/N = RT/P 
             !>  m**3/Nmole               (J/mol)*K^-1           K                      Pa         
             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )  
          ENDDO
       ENDDO


       !
       !-- Initialize
       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 )

       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 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) <= 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 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 )

             ENDIF 

          ELSE
                 
             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
        
          !
          !-- 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 /)
          
          IF ( .NOT. ALLOCATED( time_factor ) ) ALLOCATE( time_factor(1) )

          !
          !-- Get time-factor for specific hour
          index_hh = hour_of_day

          time_factor(1) = par_emis_time_factor(index_hh)

       ENDIF

       
       !
       !--  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

          !
          !-- Cycle over species:
          !-- nspec_out represents the number of species in common between the emission input data
          !-- and the chemistry mechanism used
          DO ispec=1,nspec_out  
    
             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%                                 &
                                                            preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1) * &
                                                                con_factor
          
          ENDDO
          
       !
       !-- DEFAULT case
       ELSEIF ( TRIM( mode_emis ) == "DEFAULT" )  THEN

          !
          !-- Allocate array for the emission value corresponding to a specific category and time factor
          ALLOCATE( delta_emis(nys:nyn,nxl:nxr) )  

          !
          !-- Cycle over categories
          DO icat = 1, emt_att%ncat
 
             !
             !-- Cycle over Species
             !-- nspec_out represents the number of species in common between the emission input data
             !-- and the chemistry mechanism used
             DO ispec = 1, nspec_out

                emis(nys:nyn,nxl:nxr) = emt(match_spec_input(ispec))%default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)


                !
                !-- NOx
                IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "NO" )  THEN
                
                   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
                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 
                ELSEIF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN

                   !
                   !-- Cycle over PM10 components
                   DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) )  

                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &               !<  kg/m2*s
                                                     emt_att%pm_comp(icat,i_pm_comp,3) * con_factor * hour_per_day 

                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+ &
                                                                 delta_emis(nys:nyn,nxl:nxr) 

                   ENDDO

                !
                !-- VOCs
                ELSEIF ( SIZE( match_spec_voc_input ) > 0 )  THEN

                   DO ivoc = 1, SIZE( match_spec_voc_input )

                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == TRIM( emt_att%voc_name(ivoc) ) )  THEN    

                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
                                                        emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) * &
                                                         con_factor * hour_per_day

                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emis_distribution(1,nys:nyn,nxl:nxr,ispec) + &
                                                                       delta_emis(nys:nyn,nxl:nxr) 

                      ENDIF                       

                   ENDDO
                
                !
                !-- 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
                
                emis(:,:)= 0 
                
             ENDDO
             
             delta_emis(:,:)=0 
          
          ENDDO

       ENDIF

       
!
!-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to cssws
!
!-- PARAMETERIZED mode
       !
       !-- Units of inputs are micromoles/(m**2*s)
       IF ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN

          IF ( street_type_f%from_file )  THEN

!
!--          Streets are lsm surfaces, hence, no usm surface treatment required. 
!--          However, urban surface may be initialized via default initialization
!--          in surface_mod, e.g. at horizontal urban walls that are at k == 0
!--          (building is lower than the first grid point). Hence, in order to 
!--          have only emissions at streets, set the surfaces emissions to zero
!--          at urban walls. 
             IF ( surf_usm_h%ns >=1 )  surf_usm_h%cssws = 0.0_wp
!
!--          Now, treat land-surfaces.
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
                k = surf_lsm_h%k(m)

                IF ( street_type_f%var(j,i) >= main_street_id  .AND. street_type_f%var(j,i) < max_street_id )  THEN

                   !
                   !-- Cycle over matched species
                   DO  ispec = 1, nspec_out

                      !
                      !-- PMs are already in kilograms 
                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1 "  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) )=="PM10")  THEN

                         !
                         !--           kg/(m^2*s) * kg/m^3
                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
                                                                        emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
                                                                         rho_air(k)                                    !< in kg/m^3
                         
                      !
                      !-- Other Species
                      !-- Inputs are micromoles
                      ELSE

                         !   
                         !--             ppm/s *m *kg/m^3                
                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec)) * &
                                                                        emis_distribution(1,j,i,ispec) *            &  !< in micromoles/(m^2*s)
                                                                         conv_to_ratio(k,j,i) *                     &  !< in m^3/Nmole     
                                                                          rho_air(k)                                   !< in kg/m^3

                      ENDIF

                   ENDDO


                ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND. street_type_f%var(j,i) < main_street_id )  THEN

                   !
                   !-- Cycle over matched species
                   DO  ispec = 1, nspec_out

                      !
                      !-- PMs are already in kilograms
                      IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1"  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN

                         !
                         !--           kg/(m^2*s) * kg/m^3
                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
                                                                        emis_distribution(1,j,i,ispec) *            &  !< in kg/(m^2*s)
                                                                         rho_air(k)                                    !< in kg/m^3   
                      !
                      !-- Other species
                      !-- Inputs are micromoles
                      ELSE
                    
                         !   
                         !--             ppm/s *m *kg/m^3      
                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
                                                                        emis_distribution(1,j,i,ispec) *            &  !< in micromoles/(m^2*s)
                                                                         conv_to_ratio(k,j,i) *                     &  !< in m^3/Nmole       
                                                                          rho_air(k)                                   !< in kg/m^3   
                      ENDIF

                   ENDDO

                ELSE

                   !
                   !-- If no street type is defined, then assign zero emission to all the species
                   surf_lsm_h%cssws(:,m) = 0.0_wp

                ENDIF

             ENDDO

          ENDIF

       !
       !-- For both DEFAULT and PRE-PROCESSED mode
       ELSE    
       

          DO ispec = 1, nspec_out 
                   
!
!--          Default surfaces
             DO  m = 1, surf_def_h(0)%ns

                i = surf_def_h(0)%i(m)
                j = surf_def_h(0)%j(m)

                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN

                   !
                   !-- PMs
                   IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1"  &
                        .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
                
                      !
                      !--            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 > 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) *  &       
                                                                           ( 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  
  

                      ENDIF

                   ENDIF

                ENDIF

             ENDDO
         
!
!--          LSM surfaces
             DO m = 1, surf_lsm_h%ns

                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
                k = surf_lsm_h%k(m)

                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN

                   !
                   !-- PMs
                   IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1"  &
                        .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
   
                      !
                      !--         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)            
                      surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &
                                                                     rho_air(k)                                !< in kg/m^3
 
                   ELSE

                      !
                      !-- VOCs
                      IF ( len_index_voc > 0 .AND. emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
                         !
                         !--          (ppm/s) * m * kg/m^3                        mole/(m^2/s)   
                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  &     
                                                                        conv_to_ratio(k,j,i) *           &     !< in m^3/mole
                                                                         ratio2ppm *                     &     !< in ppm
                                                                          rho_air(k)                           !< in kg/m^3  

                      !
                      !-- Other species
                      ELSE
                         !
                         !--           (ppm/s) * m  * kg/m^3                    kg/(m^2/s)
                         surf_lsm_h%cssws(match_spec_model(ispec),m) = emis_distribution(1,j,i,ispec) *  & 
                                                                        ( 1.0_wp / emt_att%xm(ispec) ) * &     !< in mole/kg
                                                                         conv_to_ratio(k,j,i) *          &     !< in m^3/mole
                                                                          ratio2ppm *                    &     !< in ppm
                                                                           rho_air(k)                          !< in kg/m^3     
                                                   
                      ENDIF

                   ENDIF

                ENDIF

             ENDDO

!
!--          USM surfaces
             DO m = 1, surf_usm_h%ns

                i = surf_usm_h%i(m)
                j = surf_usm_h%j(m)
                k = surf_usm_h%k(m)

                IF ( emis_distribution(1,j,i,ispec) > 0.0_wp )  THEN

                   !
                   !-- PMs
                   IF ( TRIM( spc_names(match_spec_model(ispec)) ) == "PM1" &
                        .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM25"  &
                          .OR. TRIM( spc_names(match_spec_model(ispec)) ) == "PM10" )  THEN
                   
                      !
                      !--          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 > 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) *  &   
                                                                        ( 1.0_wp / emt_att%xm(ispec) ) * &     !< in mole/kg
                                                                         conv_to_ratio(k,j,i) *          &     !< in m^3/mole
                                                                          ratio2ppm*                     &     !< in ppm
                                                                           rho_air(k)                          !< in kg/m^3   


                      ENDIF

                   ENDIF

                ENDIF
 
             ENDDO

          ENDDO

       ENDIF

       !
       !-- Deallocate time_factor in case of DEFAULT mode) 
       IF ( ALLOCATED ( time_factor ) ) DEALLOCATE( time_factor )

   ENDIF

 END SUBROUTINE chem_emissions_setup

 END MODULE chem_emissions_mod