Index: palm/trunk/SOURCE/basic_constants_and_equations_mod.f90
===================================================================
--- palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 4054)
+++ palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 4055)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Added rgas_univ (universal gas constant) (E.C. Chan)
+! 
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! OpenACC port for SPEC
 ! 
@@ -62,4 +66,5 @@
     REAL(wp), PARAMETER ::  pi = 3.141592654_wp                       !< PI
     !$ACC DECLARE COPYIN(pi)
+    REAL(wp), PARAMETER ::  rgas_univ = 8.31446261815324_wp           !< universal gas constant (J K-1 mol-1)
     REAL(wp), PARAMETER ::  rho_l = 1.0E3_wp                          !< density of water (kg m-3)
     REAL(wp), PARAMETER ::  rho_nacl = 2165.0_wp                      !< density of NaCl (kg m-3)
Index: palm/trunk/SOURCE/chem_emissions_mod.f90
===================================================================
--- palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 4054)
+++ palm/trunk/SOURCE/chem_emissions_mod.f90	(revision 4055)
@@ -27,4 +27,11 @@
 ! -----------------
 ! $Id$
+! - replaced local thermo. constants w/ module definitions in
+!   basic_constants_and_equations_mod (rgas_univ, p_0, r_d_cp) 
+! - initialize array emis_distribution immediately following allocation
+! - lots of minor formatting changes based on review sesson in 20190325 
+!   (E.C. Chan)
+! 
+! 3968 2019-05-13 11:04:01Z suehring
 ! - in subroutine chem_emissions_match replace all decision structures relating to
 !   mode_emis to emiss_lod
@@ -54,5 +61,5 @@
 ! 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
+! bug fix: spcs_name replaced by nvar in DO  loops
 !
 ! 3591 2018-11-30 17:37:32Z suehring
@@ -118,10 +125,13 @@
  MODULE chem_emissions_mod
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                          &
         ONLY:  rho_air
 
-    USE control_parameters,                                                    &
-        ONLY:  debug_output,                                                   &
-               end_time, message_string, initializing_actions,                 &
+    USE basic_constants_and_equations_mod,                                  &
+        ONLY:  rgas_univ, p_0, rd_d_cp
+
+    USE control_parameters,                                                 &
+        ONLY:  debug_output,                                                &
+               end_time, message_string, initializing_actions,              &
                intermediate_timestep_count, dt_3d
  
@@ -134,19 +144,19 @@
 #endif
 
-    USE netcdf_data_input_mod,                                                  &
+    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,                                        &
+    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,                                                      &
+    USE chem_gasphase_mod,                                                   &
         ONLY: nvar, spc_names
  
     USE chem_modules
 
-    USE statistics,                                                             &   
+    USE statistics,                                                          &   
         ONLY:  weight_pres
 
@@ -159,4 +169,5 @@
     CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
 
+    INTEGER(iwp) ::  dt_emis                         !< Time Step Emissions 
     INTEGER(iwp) ::  i                               !< index 1st selected dimension (some dims are not spatial)
     INTEGER(iwp) ::  j                               !< index 2nd selected dimension 
@@ -165,16 +176,11 @@
     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) ::  len_index                       !< length of index (used for several indices)
+    INTEGER(iwp) ::  len_index_pm                    !< length of PMs index
+    INTEGER(iwp) ::  len_index_voc                   !< length of voc index
     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)
+
+    REAL(wp) ::  conversion_factor                   !< Units Conversion Factor
 
     SAVE
@@ -204,5 +210,5 @@
 !
 !-- Public Variables
-    PUBLIC con_factor, len_index, len_index_pm, len_index_voc
+    PUBLIC conversion_factor, len_index, len_index_pm, len_index_voc
 
  CONTAINS
@@ -249,6 +255,4 @@
 SUBROUTINE chem_emissions_match( emt_att,len_index )   
 
-
-
     INTEGER(iwp)  ::  ind_inp                    !< Parameters for cycling through chemical input species 
     INTEGER(iwp)  ::  ind_mod                    !< Parameters for cycling through chemical model species 
@@ -282,12 +286,15 @@
           len_index = 0
 
-          IF  ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
+! number of species and number of matched species can be different
+! but call is only made if both are greater than zero
+
+          IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
 
 !
 !-- Cycle over model species
 
-             DO ind_mod = 1, nvar
+             DO  ind_mod = 1, nvar
                 ind_inp = 1
-                DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )    !< 'novalue' is the default  
+                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
@@ -310,7 +317,7 @@
                 len_index = 0
 
-                DO ind_mod = 1, nvar  
+                DO  ind_mod = 1, nvar  
                    ind_inp = 1
-                   DO WHILE  ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
+                   DO  WHILE  ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
                       IF  ( TRIM(surface_csflux_name(ind_inp)) ==          &
                             TRIM(spc_names(ind_mod))            )  THEN
@@ -326,8 +333,8 @@
 !-- 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
+                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:'             //          &
@@ -370,18 +377,22 @@
        CASE (1)
 
-          len_index = 0          ! TOTAL number of species (to be accumulated)  
-          len_index_voc = 0      ! TOTAL number of VOCs (to be accumulated)
-          len_index_pm = 3       ! TOTAL number of PMs: PM1, PM2.5, PM10.
- 
-          IF  ( nvar > 0 .AND. nspec_emis_inp > 0 )  THEN
+          len_index = 0          ! total number of species (to be accumulated)  
+          len_index_voc = 0      ! total number of VOCs (to be accumulated)
+          len_index_pm = 3       ! total number of PMs: PM1, PM2.5, PM10.
+
+!
+!-- number of model species and input species could be different
+!-- but process this only when both are non-zero
+ 
+          IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
 
 !
 !-- Cycle over model species
-             DO ind_mod = 1, nvar
+             DO  ind_mod = 1, nvar
 
 !
 !-- Cycle over input species
 
-                DO ind_inp = 1, nspec_emis_inp
+                DO  ind_inp = 1, nspec_emis_inp
 
 !
@@ -389,5 +400,5 @@
 
                    IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
-                      DO ind_voc= 1, emt_att%nvoc
+                      DO  ind_voc= 1, emt_att%nvoc
                         
                          IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
@@ -483,14 +494,14 @@
                 len_index_voc = 0
                 
-                DO ind_mod = 1, nvar
-                   DO ind_inp = 1, nspec_emis_inp 
+                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%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) ) ==                  &
@@ -617,12 +628,12 @@
           len_index_voc = 0
 
-          IF  ( nvar > 0 .AND. (nspec_emis_inp > 0) )  THEN
+          IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
 !
 !-- Cycle over model species
-             DO ind_mod = 1, nvar 
+             DO  ind_mod = 1, nvar 
 
 !
 !-- Cycle over input species  
-                DO ind_inp = 1, nspec_emis_inp
+                DO  ind_inp = 1, nspec_emis_inp
 
 !
@@ -630,5 +641,5 @@
 
                    IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN        
-                      DO ind_voc = 1, emt_att%nvoc
+                      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
@@ -674,18 +685,18 @@
 !-- Cycle over model species
 
-                DO ind_mod = 1, nvar
+                DO  ind_mod = 1, nvar
  
 !
 !-- Cycle over Input species  
 
-                   DO ind_inp = 1, nspec_emis_inp
+                   DO  ind_inp = 1, nspec_emis_inp
 
 !
 !-- VOCs
 
-                      IF  ( TRIM(emt_att%species_name(ind_inp) ) == "VOC" .AND.    &
+                      IF  ( TRIM(emt_att%species_name(ind_inp) ) == "VOC"  .AND.     &
                            ALLOCATED(match_spec_voc_input) )  THEN
                          
-                         DO ind_voc= 1, emt_att%nvoc
+                         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
@@ -716,5 +727,5 @@
 
 !
-!-- in case there are no species matching
+!-- in case there are no species matching (just informational message)
 
                 message_string = 'Non of given emission species'            //         &
@@ -767,12 +778,13 @@
 !> fluxes at timestep 0 
 !------------------------------------------------------------------------------!
+
  SUBROUTINE chem_emissions_init
 
-    USE netcdf_data_input_mod,                                                 &
+    USE netcdf_data_input_mod,                                              &
         ONLY:  chem_emis, chem_emis_att
     
     IMPLICIT NONE
  
-    INTEGER(iwp)                :: ispec                                            !< running index
+    INTEGER(iwp) :: ispec                        !< running index
 
 !    
@@ -810,5 +822,5 @@
        ALLOCATE( chem_emis_att%xm(n_matched_vars) )
 
-       DO ispec = 1, n_matched_vars
+       DO  ispec = 1, n_matched_vars
           SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
              CASE ( 'SO2'  );  chem_emis_att%xm(ispec) = xm_S + xm_O * 2
@@ -837,5 +849,4 @@
 !--        is defined, and mode if not) as we can easily take out
 !--        the case structure for mode_emis later on.
-!--        (ecc 20140424)
 
        IF   ( emiss_lod < 0 )  THEN  !-- no LOD defined (not likely)
@@ -845,5 +856,5 @@
              CASE ( 'PARAMETERIZED' )     ! LOD 0
 
-                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
                 ENDIF
@@ -853,5 +864,5 @@
              CASE ( 'DEFAULT' )           ! LOD 1
 
-                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
                 ENDIF
@@ -861,5 +872,5 @@
              CASE ( 'PRE-PROCESSED' )     ! LOD 2
 
-                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,n_matched_vars) )
                 ENDIF
@@ -875,5 +886,5 @@
              CASE ( 0 )     ! parameterized mode
 
-                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
                 ENDIF
@@ -883,5 +894,5 @@
              CASE ( 1 )     ! default mode
 
-                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,0:ny,0:nx,n_matched_vars) )
                 ENDIF
@@ -891,5 +902,5 @@
              CASE ( 2 )     ! pre-processed mode
 
-                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(nzb:nzt+1,0:ny,0:nx,n_matched_vars) )
                 ENDIF
@@ -900,4 +911,9 @@
 
        ENDIF
+
+!
+! -- initialize
+
+       emis_distribution = 0.0_wp
 
     ENDIF
@@ -917,9 +933,11 @@
  SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
  
-   USE surface_mod,                                                  &
+   USE surface_mod,                                               &
        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
-   USE netcdf_data_input_mod,                                        &
+
+   USE netcdf_data_input_mod,                                     &
        ONLY:  street_type_f
-   USE arrays_3d,                                                    &        
+
+   USE arrays_3d,                                                 &        
        ONLY: hyp, pt 
 
@@ -937,47 +955,48 @@
 
     INTEGER(iwp) ::  i                                                          !< running index for grid in x-direction
+    INTEGER(iwp) ::  i_pm_comp                                                  !< index for number of PM components
+    INTEGER(iwp) ::  icat                                                       !< Index for number of categories
+    INTEGER(iwp) ::  ispec                                                      !< index for number of species
+    INTEGER(iwp) ::  ivoc                                                       !< Index for number of VOCs
     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
-
-    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
+    REAL(wp), PARAMETER ::  hour_per_year =  8760.0_wp  !< number of hours in a year of 365 days  
+    REAL(wp), PARAMETER ::  hour_per_day  =    24.0_wp  !< number of hours in a day
+    REAL(wp), PARAMETER ::  s_per_hour    =  3600.0_wp  !< number of sec per hour (s)/(hour)    
+    REAL(wp), PARAMETER ::  s_per_day     = 86400.0_wp  !< number of sec per day (s)/(day)  
+
+    REAL(wp), PARAMETER ::  day_to_s      = 1.0_wp/s_per_day                   !< conversion day   -> sec
+    REAL(wp), PARAMETER ::  hour_to_s     = 1.0_wp/s_per_hour                  !< conversion hours -> sec
+    REAL(wp), PARAMETER ::  year_to_s     = 1.0_wp/(s_per_hour*hour_per_year)  !< conversion year  -> sec
     !
-    !-- 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: MASS
+    REAL(wp), PARAMETER ::  g_to_kg       = 1.0E-03_wp     !< Conversion from g to kg (kg/g) 
+    REAL(wp), PARAMETER ::  miug_to_kg    = 1.0E-09_wp     !< Conversion from g to kg (kg/g)
+    REAL(wp), PARAMETER ::  tons_to_kg    = 100.0_wp       !< Conversion from tons to kg (kg/tons)    
     !
-    !-- CONVERSION FACTORS: fraction to ppm
-    REAL(wp), PARAMETER ::  ratio2ppm  = 1.0e06  
-    !------------------------------------------------------    
+    !-- CONVERSION FACTORS: PPM
+    REAL(wp), PARAMETER ::  ratio2ppm     = 1.0E+06_wp 
+ 
+    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                !< temporary variable for abs. temperature
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  delta_emis   !< incremental emission factor                     
+    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  time_factor  !< factor for time scaling of emissions
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emis         !< emission factor
 
     IF  ( emission_output_required )  THEN
 
-       !
-       !-- Set emis_dt  
+!
+!-- Set emis_dt to be used - since chemistry ODEs can be stiff, the option
+!-- to solve them at every RK substep is present to help improve stability
+!-- should the need arises
+ 
        IF  ( call_chem_at_all_substeps )  THEN
 
@@ -994,24 +1013,24 @@
  !-- 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
+        IF  ( TRIM( mode_emis ) == "DEFAULT"  .OR.  TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
 
           SELECT CASE ( TRIM( emt_att%units ) )
 
-             CASE ( 'kg/m2/s',    'KG/M2/S'    );  con_factor = 1.0_wp                   ! kg
-             CASE ( 'kg/m2/hour', 'KG/M2/HOUR' );  con_factor = hour_to_s
-             CASE ( 'kg/m2/day',  'KG/M2/DAY'  );  con_factor = day_to_s
-             CASE ( 'kg/m2/year', 'KG/M2/YEAR' );  con_factor = year_to_s
-
-             CASE ( 'ton/m2/s',    'TON/M2/S'    );  con_factor = tons_to_kg             ! tonnes
-             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' );  con_factor = tons_to_kg*hour_to_s
-             CASE ( 'ton/m2/year', 'TON/M2/YEAR' );  con_factor = tons_to_kg*year_to_s
-
-             CASE ( 'g/m2/s',    'G/M2/S'    );  con_factor = g_to_kg                    ! grams
-             CASE ( 'g/m2/hour', 'G/M2/HOUR' );  con_factor = g_to_kg*hour_to_s
-             CASE ( 'g/m2/year', 'G/M2/YEAR' );  con_factor = g_to_kg*year_to_s
-
-             CASE ( 'micrograms/m2/s',    'MICROGRAMS/M2/S'    );  con_factor = miug_to_kg            ! ug
-             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' );  con_factor = miug_to_kg*hour_to_s
-             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' );  con_factor = miug_to_kg*year_to_s
+             CASE ( 'kg/m2/s', 'KG/M2/S'    );     conversion_factor = 1.0_wp                   ! kg
+             CASE ( 'kg/m2/hour', 'KG/M2/HOUR' );  conversion_factor = hour_to_s
+             CASE ( 'kg/m2/day', 'KG/M2/DAY'  );   conversion_factor = day_to_s
+             CASE ( 'kg/m2/year', 'KG/M2/YEAR' );  conversion_factor = year_to_s
+
+             CASE ( 'ton/m2/s', 'TON/M2/S'    );     conversion_factor = tons_to_kg             ! tonnes
+             CASE ( 'ton/m2/hour', 'TON/M2/HOUR' );  conversion_factor = tons_to_kg*hour_to_s
+             CASE ( 'ton/m2/year', 'TON/M2/YEAR' );  conversion_factor = tons_to_kg*year_to_s
+
+             CASE ( 'g/m2/s', 'G/M2/S'    );     conversion_factor = g_to_kg                    ! grams
+             CASE ( 'g/m2/hour', 'G/M2/HOUR' );  conversion_factor = g_to_kg*hour_to_s
+             CASE ( 'g/m2/year', 'G/M2/YEAR' );  conversion_factor = g_to_kg*year_to_s
+
+             CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S'    );     conversion_factor = miug_to_kg            ! ug
+             CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' );  conversion_factor = miug_to_kg*hour_to_s
+             CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' );  conversion_factor = miug_to_kg*year_to_s
 
 !
@@ -1031,11 +1050,12 @@
 !-- Conversion factor to convert  kg/m**2/s to ppm/s
 
-       DO i = nxl, nxr
-          DO j = nys, nyn
+       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
+             tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) *                     &
+                                       ( hyp(nzb:nzt+1) / p_0 )**rd_d_cp
  
 !            
@@ -1046,6 +1066,11 @@
 !-- V/N = RT/P 
 
+             conv_to_ratio(nzb:nzt+1,j,i) =  rgas_univ *                       &  ! J K-1 mol-1
+                                             tmp_temp(nzb:nzt+1,j,i) /         &  ! K
+                                             hyp(nzb:nzt+1)                       ! Pa
+
+! (ecc) for reference
 !                   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))) )  
+!             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )  
 
           ENDDO
@@ -1053,8 +1078,9 @@
 
 
+! (ecc) moved initialization immediately after allocation
 !
 !-- Initialize
 
-       emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
+!       emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
 
   
@@ -1084,5 +1110,5 @@
 !-- Allocate array where to store temporary emission values
 
-          IF  ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
+          IF  (  .NOT.  ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
 
 !
@@ -1129,5 +1155,5 @@
 !
 !-- Update time indices
-             CALL time_default_indices( daytype_mdh, month_of_year, day_of_month,     &
+             CALL time_default_indices( daytype_mdh, month_of_year, day_of_month,  &
                   hour_of_day, index_mm, index_dd,index_hh )
 
@@ -1178,10 +1204,10 @@
 !-- assign constant values of time factors, diurnal profile for traffic sector
 
-          par_emis_time_factor( : ) =  (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039,      &
-                                          0.056, 0.053, 0.051, 0.051, 0.052, 0.055,      &
-                                          0.059, 0.061, 0.064, 0.067, 0.069, 0.069,      &
+          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))
+          IF  (  .NOT.  ALLOCATED (time_factor) )  ALLOCATE (time_factor(1))
 
 !
@@ -1205,5 +1231,5 @@
 !       IF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
 
-          DO ispec = 1, n_matched_vars
+          DO  ispec = 1, n_matched_vars
 
 !
@@ -1214,5 +1240,5 @@
                        time_factor(1) * hour_to_s
 
-          ENDDO 
+          ENDDO  
 
           
@@ -1234,5 +1260,5 @@
 !-- Cycle over categories
 
-          DO icat = 1, emt_att%ncat
+          DO  icat = 1, emt_att%ncat
  
 !
@@ -1240,5 +1266,5 @@
 !-- in common between the emission input data and the chemistry mechanism used
 
-             DO ispec = 1, n_matched_vars
+             DO  ispec = 1, n_matched_vars
 
                 emis(nys:nyn,nxl:nxr) =                    &
@@ -1254,5 +1280,5 @@
                                                  time_factor(icat) *           &
                                                  emt_att%nox_comp(icat,1) *    &
-                                                 con_factor * hour_per_day
+                                                 conversion_factor * hour_per_day
 
                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
@@ -1267,5 +1293,5 @@
                                                  time_factor(icat) *           &
                                                  emt_att%nox_comp(icat,2) *    &
-                                                 con_factor * hour_per_day
+                                                 conversion_factor * hour_per_day
 
                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
@@ -1280,5 +1306,5 @@
                                                  time_factor(icat) *           &
                                                  emt_att%sox_comp(icat,1) *    &
-                                                 con_factor * hour_per_day
+                                                 conversion_factor * hour_per_day
 
                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
@@ -1295,5 +1321,5 @@
                                                  time_factor(icat) *           &
                                                  emt_att%sox_comp(icat,2) *    &
-                                                 con_factor * hour_per_day
+                                                 conversion_factor * hour_per_day
 
                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
@@ -1312,5 +1338,5 @@
                                                     time_factor(icat) *                  &
                                                     emt_att%pm_comp(icat,i_pm_comp,1) *  &
-                                                    con_factor * hour_per_day 
+                                                    conversion_factor * hour_per_day 
 
                       emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
@@ -1330,5 +1356,5 @@
                                                     time_factor(icat) *                  &
                                                     emt_att%pm_comp(icat,i_pm_comp,2) *  &
-                                                    con_factor * hour_per_day 
+                                                    conversion_factor * hour_per_day 
 
                       emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
@@ -1348,5 +1374,5 @@
                                                     time_factor(icat)     *              &
                                                     emt_att%pm_comp(icat,i_pm_comp,3) *  &
-                                                    con_factor * hour_per_day 
+                                                    conversion_factor * hour_per_day 
 
                       emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                       &
@@ -1369,5 +1395,5 @@
                                                        time_factor(icat) *               &
                                                        emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) *   &
-                                                       con_factor * hour_per_day
+                                                       conversion_factor * hour_per_day
 
                          emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                    &
@@ -1386,5 +1412,5 @@
                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *                 &
                                                  time_factor(icat) *                     &
-                                                 con_factor * hour_per_day
+                                                 conversion_factor * hour_per_day
  
                    emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                          &
@@ -1394,9 +1420,9 @@
                 ENDIF  ! TRIM spc_names
                 
-                emis(:,:)= 0 
+                emis = 0 
                 
              ENDDO
              
-             delta_emis(:,:)=0 
+             delta_emis = 0 
           
           ENDDO
@@ -1414,5 +1440,5 @@
 !-- in common between the emission input data and the chemistry mechanism used
 
-          DO ispec = 1, n_matched_vars  
+          DO  ispec = 1, n_matched_vars  
  
 ! (ecc)   
@@ -1420,5 +1446,5 @@
                        emt(match_spec_input(ispec))%                                     &
                            preproc_emission_data(index_hh,1,nys+1:nyn+1,nxl+1:nxr+1) *   &
-                       con_factor
+                       conversion_factor
 
 
@@ -1426,5 +1452,5 @@
 !                       emt(match_spec_input(ispec))%                                     &
 !                           preproc_emission_data(index_hh,1,:,:) *   &
-!                       con_factor
+!                       conversion_factor
           ENDDO
 
@@ -1464,4 +1490,9 @@
                 k = surf_lsm_h%k(m)
 
+!
+!-- set everything to zero then reassign according to street type
+
+                surf_lsm_h%cssws(:,m) = 0.0_wp
+
                 IF  ( street_type_f%var(j,i) >= main_street_id    .AND.                  &
                       street_type_f%var(j,i) < max_street_id   )  THEN
@@ -1475,6 +1506,6 @@
 !-- PMs are already in kilograms 
 
-                      IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"     .OR. &
-                            TRIM(spc_names(match_spec_model(ispec))) == "PM25"    .OR. &
+                      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
 
@@ -1516,6 +1547,6 @@
 !-- PMs are already in kilograms
 
-                      IF  ( TRIM(spc_names(match_spec_model(ispec))) == "PM1"     .OR. &
-                            TRIM(spc_names(match_spec_model(ispec))) == "PM25"    .OR. &
+                      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
 
@@ -1548,7 +1579,8 @@
 !-- If no street type is defined, then assign zero emission to all the species
 
-                ELSE
-
-                   surf_lsm_h%cssws(:,m) = 0.0_wp
+! (ecc) moved to front (for reference)
+!                ELSE
+!
+!                   surf_lsm_h%cssws(:,m) = 0.0_wp
 
                 ENDIF  ! street type
@@ -1566,5 +1598,5 @@
        
 
-          DO ispec = 1, n_matched_vars 
+          DO  ispec = 1, n_matched_vars 
                    
 !
@@ -1702,5 +1734,5 @@
 !
 !-- VOCs
-                      IF  ( len_index_voc > 0                                           .AND.      &
+                      IF  ( len_index_voc > 0                                         .AND.      &
                             emt_att%species_name(match_spec_input(ispec)) == "VOC" )  THEN
 
