Index: palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3373)
+++ palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 3458)
@@ -27,4 +27,15 @@
 ! -----------------
 ! $Id$
+! 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
 ! 
@@ -40,7 +51,7 @@
 ! Authors:
 ! --------
-! Emmanuele Russo Fu-Berlin
-! Sabine Banzhaf  FU-Berlin
-! Martijn Schaap  FU-Berlin, TNO Utrecht
+! @author Emmanuele Russo (Fu-Berlin)
+! @author Sabine Banzhaf  (FU-Berlin)
+! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
 ! 
 ! Description:
@@ -49,6 +60,7 @@
 !>
 !> @todo Check_parameters routine should be developed: for now it includes just one condition
-!> @todo Add option for capital or small letters in the matching routine
 !> @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
 !> @note <Enter notes on the module>
 !> @bug  <Enter known bugs here>
@@ -61,5 +73,5 @@
 
     USE control_parameters,                                                    &
-       ONLY:  initializing_actions,end_time, message_string,                   &
+       ONLY:  initializing_actions, end_time, message_string,                  &
               intermediate_timestep_count, dt_3d
  
@@ -67,14 +79,17 @@
 
     USE kinds
+
+#if defined ( __netcdf )
 
     USE netcdf_data_input_mod,                                                  &
        ONLY: chem_emis_att_type, chem_emis_val_type
 
-#if defined ( __netcdf )
     USE NETCDF
+
 #endif
 
     USE date_and_time_mod,                                                      &
-       ONLY: time_default_indices, month_of_year, day_of_month, hour_of_day,    &
+       ONLY: time_default_indices, time_preprocessed_indices,                  &
+             month_of_year, day_of_month, hour_of_day,                          &
              index_mm, index_dd, index_hh
 
@@ -92,7 +107,5 @@
     
     CHARACTER (LEN=80)                               :: filename_emis                   !> Variable for the name of the netcdf input file
-    CHARACTER (LEN=80), ALLOCATABLE, DIMENSION(:)    :: string_values                   !> Output of string variables read from netcdf input file
-
-    INTEGER(iwp)                                     :: n_dims                          !> Number of dimensions of input variable netcdf file
+
     INTEGER(iwp)                                     :: i                               !> index 1st selected dimension (some dims are not spatial)
     INTEGER(iwp)                                     :: j                               !> index 2nd selected dimension 
@@ -103,9 +116,4 @@
     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), ALLOCATABLE, DIMENSION(:)          :: id_dims                         !> id of dimension of selected variable netcdf input
-    INTEGER(iwp), ALLOCATABLE, DIMENSION(:)          :: dum_var                         !> value of variable read from netcdf input
-    INTEGER(iwp)                                     :: id_ncfile                       !> id netcdf file 
-    INTEGER(iwp)                                     :: errno                           !> error number NF90_???? function
-    INTEGER(iwp)                                     :: id_var                          !< variable id
     INTEGER(iwp)                                     :: dt_emis                         !> Time Step Emissions 
     INTEGER(iwp)                                     :: len_index                       !> length of index (used for several indices)
@@ -124,38 +132,4 @@
     !> Dobson units:
     REAL, PARAMETER        ::  Dobs = 2.68668e16    ! (mlc/cm2) / DU
-   
-    ! molar weights of components
-  
-    !> atomic 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
-    REAL, PARAMETER        ::  xm_Rh    =  102.90550e-3               ! kg/mol
-  
-    !> molecule 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
-    
-    !> mass of air
-    REAL, PARAMETER        ::  xm_air   =  28.964e-3                  ! kg/mol
-    REAL, PARAMETER        ::  xmair    =  28.94                      ! g/mol; old name!
 
     !> sesalt composition:
@@ -263,6 +237,6 @@
 !-- Public Variables
 
-    PUBLIC id_ncfile,errno,id_dims,n_dims,dum_var, con_factor,   &
-           len_index,len_index_voc,len_index_pm, string_values
+    PUBLIC con_factor, len_index,len_index_voc,len_index_pm
+
  CONTAINS
 
@@ -275,5 +249,5 @@
 
 
-!TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?
+    !TBD: Where Should we put the call to chem_emissions_check_parameters? In chem_init or in check_parameters?
 
     IMPLICIT NONE
@@ -548,5 +522,6 @@
                 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(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
                 ENDIF
 
@@ -661,5 +636,5 @@
                                  ' DO NOT MATCH'                                //          &
                                  ' model chemical species'                      //          &
-                                 ' Chemistry Emissions routine is not called'         !TBD: IMPORTANT:Add a condition in chem init that sets the emission_output_required to false when len_index==0 so that the chem_emissions modules are not used
+                                 ' Chemistry Emissions routine is not called'         
                 CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 
 
@@ -724,8 +699,7 @@
                 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:'            //           &
+                   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:'            //          &
                                        ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
                                        ' but values of scaling factors for street types'    //          &
@@ -780,5 +754,5 @@
  SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out)
 
-
+#if defined( __netcdf )
 
     USE surface_mod,                                                           &
@@ -799,5 +773,5 @@
  
     INTEGER(iwp)                                                      :: ispec     !> Index to go through the emission chemical species
-#if defined( __netcdf )
+
 
 !-- Actions for initial runs : TBD: needs to be updated
@@ -854,4 +828,15 @@
      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")
@@ -864,20 +849,10 @@
            CALL chem_emissions_setup(emt_att,emt,nspec_out)
 
-
-        !> 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)
-
+        !> 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 EXP all. 
+           ! 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))
 
@@ -919,5 +894,5 @@
  IMPLICIT NONE
 
-
+#if defined( __netcdf )
   
     !--- IN/OUT
@@ -949,4 +924,7 @@
     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
@@ -964,12 +942,19 @@
     ! --- 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?
+    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 = 0.0002777778  !  (hour)/(s)
+    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.157407e-05   !  (day)/(s)
+    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 = 3.170979e-08  !  (year)/(s)
+    REAL, PARAMETER   ::  year_to_s = 1/(s_per_hour*hour_per_year)  !  (year)/(s)
 
     !-CONVERSION FACTORS: WEIGHT
@@ -984,8 +969,8 @@
     REAL(wp), PARAMETER   ::  ratio2ppm  = 1.0e06_wp  
     !------------------------------------------------------    
-#if defined( __netcdf )
+
     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 calle: for now one hour. It should be made changeable.
+    !>  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 ( call_chem_at_all_substeps )  THEN
@@ -1093,5 +1078,7 @@
  
           message_string = 'No Units conversion required for units of chemistry emissions' // &
-                           ' of the PARAMETERIZED mode: units have to be provided in mole/m**2/s '
+                           ' 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 ) 
 
@@ -1110,5 +1097,5 @@
              !  V/N=RT/P 
 
-             !>    m**3/Nmole             (J/mol)*K^-1           K                      Pa         
+             !>    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
@@ -1125,4 +1112,7 @@
     !> 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. )
@@ -1194,4 +1184,8 @@
                                  ' 24 values for every day of the year ', .FALSE. )
         
+          !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))
@@ -1200,5 +1194,5 @@
           index_hh=hour_of_day
 
-          time_factor(1)=emt_att%par_emis_time_factor(index_hh)
+          time_factor(1) = par_emis_time_factor(index_hh)
 
        ENDIF
@@ -1212,6 +1206,6 @@
           DO ispec=1,nspec_out
 
-             !> Values are still micromoles/(m**2*s). Units are not given by the user for this mode and are always micromoles (or micrograms for PMs)
-             emis_distribution(1,nys:nyn,nxl:nxr,ispec)=surface_csflux(match_spec_input(ispec))*time_factor(1)
+             !> 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
 
           ENDDO 
@@ -1223,9 +1217,9 @@
           DO ispec=1,nspec_out !> nspec_out represents the number of species in common between 
                                !  the emission input data and the chemistry mechanism used
-
-             emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))%                                  &
-                                                                      preproc_emission_data(index_hh,nzb:nzt+1,nys:nyn,nxl:nxr)* &
-                                                                   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
 
@@ -1249,38 +1243,34 @@
 !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                       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
+                !>             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)
+                   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
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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
-                      
-                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
-                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   !>             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
-   
-                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
-                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   !>             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)
  
 
@@ -1292,10 +1282,8 @@
                    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
+                      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)+delta_emis(nys:nyn,nxl:nxr)
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
  
                    ENDDO
@@ -1307,10 +1295,8 @@
                    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
+                      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)
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
  
                    ENDDO
@@ -1322,10 +1308,8 @@
                    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
+                      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 
                                                                                                  
 
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
-                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 
 
                    ENDDO
@@ -1339,9 +1323,8 @@
                       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
-
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                         &
-                               emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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                       
@@ -1352,8 +1335,7 @@
                 ELSE
 
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor
- 
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                               &
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*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==)
@@ -1396,7 +1378,6 @@
                       DO  ispec=1,nspec_out
 
-                         !> PMs are already in mass units:micrograms:they have to be converted to kilograms 
-                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"         &
-                             .OR.  TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
+                         !> 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
 
@@ -1412,9 +1393,9 @@
                          !> 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))*  &
+                            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)*                        &
+                                                                          emis_distribution(1,j,i,ispec) *              &
+                         !                                                    m^3/Nmole
+                                                                          conv_to_ratio(k,j,i)*                         &       
                          !                                                    kg/m^3
                                                                           rho_air(k)   
@@ -1426,6 +1407,5 @@
 
 
-                   ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.                                            &
-                            street_type_f%var(j,i) < main_street_id )  THEN
+                   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
@@ -1433,10 +1413,9 @@
 
                          !> 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"  &
+                         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)) *  &
+                            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)*        &
@@ -1448,9 +1427,9 @@
                          !>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)) * &
+                            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)*                        &
+                                                                          emis_distribution(1,j,i,ispec) *              &
+                         !                                                    m^3/Nmole
+                                                                          conv_to_ratio(k,j,i)*                         &       
                          !                                                    kg/m^3
                                                                           rho_air(k)   
@@ -1486,43 +1465,47 @@
                    j = surf_def_h(0)%j(m)
 
-                   !> 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
+                   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
                    
-                      !            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(k)   
- 
- 
-                   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(k,j,i)*ratio2ppm *      &
-                      !                                                    kg/m^3
-                                                                          rho_air(k)   
-
-
-                      !> OTHER SPECIES
+                         !            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)   
+ 
+ 
                       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(k,j,i)*ratio2ppm*       &
-                         !                                                  kg/m^3
-                                                                          rho_air(k)   
-  
+                         !> 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
+                         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)   
+  
+
+                         ENDIF
 
                       ENDIF
@@ -1543,45 +1526,48 @@
                    k = surf_lsm_h%k(m)
 
-                   !> 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
-
-                      !         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)   
- 
-                   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
+                   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
+   
+                         !         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)   
+ 
                       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)   
+                         !> 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
+                         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)   
                                                       
+                         ENDIF
+
                       ENDIF
 
                    ENDIF
+
                 ENDDO
 
@@ -1598,40 +1584,42 @@
                    k = surf_usm_h%k(m)
 
-
-                   !> 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)   
-
- 
-                   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
+                   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)*        & 
+                         !                                              kg/m^3
                                                                        rho_air(k)   
 
-                      !> 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)   
-
+                         !> 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
+                         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)   
+
+
+                         ENDIF
 
                       ENDIF
