Changeset 3458 for palm/trunk

Timestamp:

Oct 30, 2018 2:51:23 PM (6 years ago)

Author:

kanani

Message:

Reintegrated fixes/changes from branch chemistry

Location:

palm/trunk

Files:

• SOURCE/Makefile (modified) (2 diffs)
• SOURCE/chem_emissions_mod.f90 (modified) (42 diffs)
• SOURCE/chem_modules.f90 (modified) (4 diffs)
• SOURCE/chemistry_model_mod.f90 (modified) (14 diffs)
• SOURCE/date_and_time_mod.f90 (modified) (11 diffs)
• SOURCE/flow_statistics.f90 (modified) (2 diffs)
• SOURCE/init_3d_model.f90 (modified) (2 diffs)
• SOURCE/netcdf_data_input_mod.f90 (modified) (4 diffs)
• SOURCE/palm.f90 (modified) (3 diffs)
• SOURCE/prognostic_equations.f90 (modified) (3 diffs)
• UTIL/chemistry/gasphase_preproc/kpp4palm/Makefile (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh (modified) (7 diffs)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C (modified) (3 diffs)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.h (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.h (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.C (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.h (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file_vec.C (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/main.C (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.C (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.h (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/utils.C (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/src/utils.h (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/kpp4palm/templates/kp4_compress_header (added)
• UTIL/chemistry/gasphase_preproc/kpp4palm/templates/kp4_compress_subroutines (added)
• UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_cbm4/chem_gasphase_mod.f90 (modified) (12 diffs)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp (added)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.f90 (added)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.kpp (added)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/phstatp.eqn (added)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/phstatp.spc (added)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90 (modified) (2 diffs)
• UTIL/chemistry/gasphase_preproc/run_kpp4palm.ksh (modified) (1 diff)

palm/trunk/SOURCE/Makefile

@@ -25 +25 @@
 # -----------------
 # $Id$
+# from chemistry branch r3443, banzhafs, Russo, forkel, basit:
+# radiation_model_mod added to chemistry_model_mod
+# Added some missing dependencies and replaced blanks with tabs
+# Added chemistry emission module
+# chemistry_model_mod added to flow_statistics
+# 
+# 3448 2018-10-29 18:14:31Z kanani
 # Adjustment of biometeorology dependencies
 # 
@@ -856 +863 @@
         modules.o \
         netcdf_data_input_mod.o \
+        radiation_model_mod.o \
         surface_mod.o \
         user_module.o

palm/trunk/SOURCE/chem_emissions_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $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 +51 @@
 ! 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 +60 @@
 !>
 !> @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 +73 @@
 
     USE control_parameters,                                                    &
-       ONLY:  initializing_actions,end_time, message_string,                   &
+       ONLY:  initializing_actions, end_time, message_string,                  &
               intermediate_timestep_count, dt_3d
  
@@ -67 +79 @@
 
     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 +107 @@
     
     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 +116 @@
     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 +132 @@
     !> 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 +237 @@
 !-- 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 +249 @@
 
 
-!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 +522 @@
                 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 +636 @@
                                  ' 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 +699 @@
                 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 +754 @@
  SUBROUTINE chem_emissions_init(emt_att,emt,nspec_out)
 
-
+#if defined( __netcdf )
 
     USE surface_mod,                                                           &
@@ -799 +773 @@
  
     INTEGER(iwp)                                                      :: ispec     !> Index to go through the emission chemical species
-#if defined( __netcdf )
+
 
 !-- Actions for initial runs : TBD: needs to be updated
@@ -854 +828 @@
      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 +849 @@
            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 +894 @@
  IMPLICIT NONE
 
-
+#if defined( __netcdf )
   
     !--- IN/OUT
@@ -949 +924 @@
     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 +942 @@
     ! --- 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 +969 @@
     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 +1078 @@
  
           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 +1097 @@
              !  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 +1112 @@
     !> 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 +1184 @@
                                  ' 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 +1194 @@
           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 +1206 @@
           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 +1217 @@
           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 +1243 @@
 !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 +1282 @@
                    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 +1295 @@
                    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 +1308 @@
                    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 +1323 @@
                       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 +1335 @@
                 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 +1378 @@
                       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 +1393 @@
                          !> 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 +1407 @@
 
 
-                   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 +1413 @@
 
                          !> 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 +1427 @@
                          !>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 +1465 @@
                    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 +1526 @@
                    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 +1584 @@
                    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

palm/trunk/SOURCE/chem_modules.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443:
+! ??
+! 
+! 3298 2018-10-02 12:21:11Z kanani
 ! - Minor formatting
 ! - Introduced Variables for Chemistry emissions Module (Russo)
@@ -43 +47 @@
 ! @author Farah Kanani-Suehring
 ! @author Basit Khan
+! @author Sabine Banzhaf
+! @author Emmanuele Russo
 !
 !------------------------------------------------------------------------------!
@@ -88 +94 @@
     LOGICAL       :: do_emis                               = .FALSE.               !< Flag for turning on chemistry emissions 
     LOGICAL       :: cs_pr_namelist_found                  = .FALSE.               !< Namelist parameter: Names of t
+    LOGICAL       :: do_depo                               = .FALSE.               !< namelist parameter for activation of deposition calculation
 
 
@@ -161 +168 @@
                                                                                    ! matching between the model and the input files
 
+    
+    !-- Selected atomic/molecular weights:
+    
+    REAL, PARAMETER        ::  xm_H     =    1.00790e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_N     =   14.00670e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_C     =   12.01115e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_S     =   32.06400e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_O     =   15.99940e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_F     =   18.99840e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Na    =   22.98977e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Cl    =   35.45300e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Rn222 =  222.00000e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Pb210 =  210.00000e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Ca    =   40.07800e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_K     =   39.09800e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Mg    =   24.30500e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Pb    =  207.20000e-3           ! kg/mol
+    REAL, PARAMETER        ::  xm_Cd    =  112.41000e-3           ! kg/mol
+    
+    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
+        
+    ! dummy weight, used for complex molecules:
+    REAL, PARAMETER        ::  xm_dummy =  1000.0e-3              ! kg/mol
+
+    
     SAVE
  END MODULE chem_modules

palm/trunk/SOURCE/chemistry_model_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443, banzhafs, basit:
+! replace surf_lsm_h%qv1(m) by q(k,j,i) for mixing ratio in chem_depo
+! bug fix in chem_depo: allow different surface fractions for one
+! surface element and set lai to zero for non vegetated surfaces
+! bug fixed in chem_data_output_2d
+! bug fix in chem_depo subroutine
+! added code on deposition of gases and particles
+! removed cs_profile_name from chem_parin
+! bug fixed in output profiles and code cleaned
+! 
+! 3449 2018-10-29 19:36:56Z suehring
 ! additional output - merged from branch resler
 ! 
@@ -151 +162 @@
 ! @author Klaus Ketelsen
 ! @author Basit Khan
+! @author Sabine Banzhaf
 !
 !
@@ -197 +209 @@
                                  initializing_actions, message_string,                             &
                                  omega, tsc, intermediate_timestep_count,                          &
-                                 intermediate_timestep_count_max,                                  &
+                                 intermediate_timestep_count_max, max_pr_user,                     &
                                  timestep_scheme, use_prescribed_profile_data, ws_scheme_sca  
    USE arrays_3d,          ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
@@ -271 +283 @@
                                          ! 'fastj'  (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
 
+    
+    !-----------------------------------------------------------------------
+    ! Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
+    !
+    INTEGER(iwp), PARAMETER :: nlu_dep = 15                   ! Number of DEPAC landuse classes (lu's)
+    INTEGER(iwp), PARAMETER :: ncmp = 10                      ! Number of DEPAC gas components
+    !------------------------------------------------------------------------
+    ! DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                              
+    !
+    INTEGER(iwp)  ::  ilu_grass              = 1        
+    INTEGER(iwp)  ::  ilu_arable             = 2        
+    INTEGER(iwp)  ::  ilu_permanent_crops    = 3         
+    INTEGER(iwp)  ::  ilu_coniferous_forest  = 4         
+    INTEGER(iwp)  ::  ilu_deciduous_forest   = 5         
+    INTEGER(iwp)  ::  ilu_water_sea          = 6        
+    INTEGER(iwp)  ::  ilu_urban              = 7        
+    INTEGER(iwp)  ::  ilu_other              = 8  
+    INTEGER(iwp)  ::  ilu_desert             = 9  
+    INTEGER(iwp)  ::  ilu_ice                = 10 
+    INTEGER(iwp)  ::  ilu_savanna            = 11 
+    INTEGER(iwp)  ::  ilu_tropical_forest    = 12 
+    INTEGER(iwp)  ::  ilu_water_inland       = 13 
+    INTEGER(iwp)  ::  ilu_mediterrean_scrub  = 14 
+    INTEGER(iwp)  ::  ilu_semi_natural_veg   = 15 
+    ! 
+    !--------------------------------------------------------------------------
+    ! NH3/SO2 ratio regimes:
+    INTEGER, PARAMETER  ::  iratns_low  = 1       ! low ratio NH3/SO2
+    INTEGER, PARAMETER  ::  iratns_high = 2       ! high ratio NH3/SO2
+    INTEGER, PARAMETER  ::  iratns_very_low = 3   ! very low ratio NH3/SO2
+    ! Default:
+    INTEGER, PARAMETER  ::  iratns_default = iratns_low
+    !
+    ! Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) naar W m-2
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: alpha   =(/0.009,0.009, 0.009,0.006,0.006, -999., -999.,0.009,-999.,-999.,0.009,0.006,-999.,0.009,0.008/)*4.57
+    !
+    ! Set temperatures per land use for F_temp
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmin =   (/12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999., 12.0,  0.0, -999., 12.0,  8.0/)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Topt =   (/26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999., 26.0, 20.0, -999., 26.0, 24.0 /)
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmax =   (/40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999., 40.0, 35.0, -999., 40.0, 39.0 /)
+    !
+    ! Set F_min:
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: F_min   =(/0.01, 0.01,  0.01, 0.1,  0.1,   -999., -999.,0.01, -999.,-999.,0.01,0.1,-999.,0.01, 0.04/)
+    
+    ! Set maximal conductance (m/s)
+    ! (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
+    ! Could be refined to a function of T and P. in Jones
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max   =(/270., 300.,  300., 140., 150.,  -999., -999.,270., -999.,-999.,270., 150.,-999.,300., 422./)/41000
+    !
+    ! Set max, min for vapour pressure deficit vpd;
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_max =(/1.3,  0.9,   0.9,  0.5,  1.0,   -999., -999.,1.3,  -999.,-999.,1.3,1.0,  -999.,0.9, 2.8/) 
+    REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_min =(/3.0,  2.8,   2.8,  3.0,  3.25,  -999., -999.,3.0,  -999.,-999.,3.0,3.25, -999.,2.8, 4.5/) 
+    !
+    !  
+    !------------------------------------------------------------------------
+
+    
     PUBLIC nest_chemistry
     PUBLIC nspec
@@ -357 +426 @@
     END INTERFACE chem_wrd_local
 
+    INTERFACE chem_depo
+       MODULE PROCEDURE chem_depo 
+    END INTERFACE chem_depo
+
+    INTERFACE drydepos_gas_depac
+       MODULE PROCEDURE drydepos_gas_depac 
+    END INTERFACE drydepos_gas_depac
+
+    INTERFACE rc_special
+       MODULE PROCEDURE rc_special 
+    END INTERFACE rc_special
+
+    INTERFACE  rc_gw
+       MODULE PROCEDURE rc_gw 
+    END INTERFACE rc_gw
+
+    INTERFACE rw_so2 
+       MODULE PROCEDURE rw_so2  
+    END INTERFACE rw_so2
+
+    INTERFACE rw_nh3_sutton
+       MODULE PROCEDURE rw_nh3_sutton 
+    END INTERFACE rw_nh3_sutton
+
+    INTERFACE rw_constant
+       MODULE PROCEDURE rw_constant 
+    END INTERFACE rw_constant
+
+    INTERFACE rc_gstom
+       MODULE PROCEDURE rc_gstom 
+    END INTERFACE rc_gstom
+
+    INTERFACE rc_gstom_emb
+       MODULE PROCEDURE rc_gstom_emb 
+    END INTERFACE rc_gstom_emb
+
+    INTERFACE par_dir_diff
+       MODULE PROCEDURE par_dir_diff 
+    END INTERFACE par_dir_diff
+
+    INTERFACE rc_get_vpd
+       MODULE PROCEDURE rc_get_vpd 
+    END INTERFACE rc_get_vpd
+
+    INTERFACE rc_gsoil_eff
+       MODULE PROCEDURE rc_gsoil_eff 
+    END INTERFACE rc_gsoil_eff
+
+    INTERFACE rc_rinc
+       MODULE PROCEDURE rc_rinc 
+    END INTERFACE rc_rinc
+
+    INTERFACE rc_rctot
+       MODULE PROCEDURE rc_rctot  
+    END INTERFACE rc_rctot
+
+    INTERFACE rc_comp_point_rc_eff
+       MODULE PROCEDURE rc_comp_point_rc_eff 
+    END INTERFACE rc_comp_point_rc_eff
+
+    INTERFACE drydepo_aero_zhang_vd
+       MODULE PROCEDURE drydepo_aero_zhang_vd 
+    END INTERFACE drydepo_aero_zhang_vd
+
+    INTERFACE get_rb_cell
+       MODULE PROCEDURE  get_rb_cell
+    END INTERFACE get_rb_cell
+
+    
 
     PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
@@ -364 +502 @@
            chem_init_profiles, chem_integrate, chem_parin,                      &
            chem_prognostic_equations, chem_rrd_local,                           &
-           chem_statistics, chem_swap_timelevel, chem_wrd_local 
+           chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo 
 
 
@@ -838 +976 @@
 !            It is possible  to plant PM10 and PM25 into the gasphase chemistry code 
 !            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
-!            set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
+!            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
              IF (spec_name(1:2) == 'PM')   THEN  
                unit = 'kg m-3'
@@ -852 +990 @@
 
 
-      RETURN
+       RETURN
     END SUBROUTINE chem_check_data_output
 !
@@ -1072 +1210 @@
                          DO  k = nzb_do, nzt_do
                               local_pf(i,j,k) = MERGE(                                                &
-                                                 chem_species(lsp)%conc(k,j,i),                       &
+                                                 chem_species(lsp)%conc_av(k,j,i),                    &
                                                  REAL( fill_value, KIND = wp ),                       &
                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
@@ -1104 +1242 @@
        USE surface_mod
 
-       !USE chem_modules, ONLY: nvar
 
        IMPLICIT NONE
@@ -1239 +1376 @@
 
        spec_name = TRIM( variable(4:) )
-       !av == 0
 
        DO lsp=1,nspec
@@ -1870 +2006 @@
                                         cs_name,                          &
                                         cs_profile,                       &
-                                        cs_profile_name,                  & 
                                         cs_surface,                       &
                                         decycle_chem_lr,                  &
                                         decycle_chem_ns,                  &            
-                                        decycle_method,                   & 
+                                        decycle_method,                   &
+                                        do_depo,                          &
                                         emiss_factor_main,                &
                                         emiss_factor_side,                &                      
@@ -2365 +2501 @@
                 DO lpr = 1, cs_pr_count
 
-                   sums_l(k,pr_palm+lpr,tn) = sums_l(k,pr_palm+lpr,tn) +    &
+                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
                                                  chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
                                                  rmask(j,i,sr)  *                                   &
@@ -2432 +2568 @@
 
     END SUBROUTINE chem_wrd_local
+
+
+
+!-------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine to calculate the deposition of gases and PMs. For now deposition
+!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
+!> considered. The deposition of particles is derived following Zhang et al.,
+!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
+!>     
+!> @TODO: Consider deposition on vertical surfaces    
+!> @TODO: Consider overlaying horizontal surfaces
+!> @TODO: Check error messages
+!-------------------------------------------------------------------------------!
+     
+    SUBROUTINE chem_depo (i,j)
+
+
+      USE control_parameters,                                                 &   
+           ONLY:  dt_3d, intermediate_timestep_count, latitude
+
+      USE arrays_3d,                                                          &
+           ONLY:  dzw, rho_air_zw
+
+      USE date_and_time_mod,                                                  &
+           ONLY: day_of_year
+
+      USE surface_mod,                                                          &
+           ONLY: surf_lsm_h, surf_usm_h, surf_type, ind_veg_wall,               &
+           ind_pav_green, ind_wat_win
+
+      USE radiation_model_mod,                                                   &
+           ONLY: zenith
+
+
+
+      IMPLICIT NONE
+
+      INTEGER,INTENT(IN)       :: i,j
+
+      INTEGER(iwp) ::  k                                        ! matching k to surface m at i,j
+      INTEGER(iwp) ::  lsp                                      ! running index for chem spcs.
+      INTEGER(iwp) ::  lu                                       ! running index for landuse classes
+      INTEGER(iwp) ::  lu_palm                                  ! index of PALM LSM vegetation_type at current surface element
+      INTEGER(iwp) ::  lup_palm                                 ! index of PALM LSM pavement_type at current surface element
+      INTEGER(iwp) ::  luw_palm                                 ! index of PALM LSM water_type at current surface element
+      INTEGER(iwp) ::  luu_palm                                 ! index of PALM USM walls/roofs at current surface element
+      INTEGER(iwp) ::  lug_palm                                 ! index of PALM USM green walls/roofs at current surface element
+      INTEGER(iwp) ::  lud_palm                                 ! index of PALM USM windows at current surface element
+      INTEGER(iwp) ::  lu_dep                                   ! matching DEPAC LU to lu_palm
+      INTEGER(iwp) ::  lup_dep                                  ! matching DEPAC LU to lup_palm
+      INTEGER(iwp) ::  luw_dep                                  ! matching DEPAC LU to luw_palm
+      INTEGER(iwp) ::  luu_dep                                  ! matching DEPAC LU to luu_palm
+      INTEGER(iwp) ::  lug_dep                                  ! matching DEPAC LU to lug_palm
+      INTEGER(iwp) ::  lud_dep                                  ! matching DEPAC LU to lud_palm
+      INTEGER(iwp) ::  m                                        ! index for horizontal surfaces
+
+      INTEGER(iwp) ::  pspec                                    ! running index
+      INTEGER(iwp) ::  i_pspec                                  ! index for matching depac gas component
+
+
+      ! Vegetation                                              !Match to DEPAC
+      INTEGER(iwp)  ::  ind_lu_user = 0                         ! --> ERROR  
+      INTEGER(iwp)  ::  ind_lu_b_soil = 1                       ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lu_mixed_crops = 2                  ! --> ilu_arable
+      INTEGER(iwp)  ::  ind_lu_s_grass = 3                      ! --> ilu_grass
+      INTEGER(iwp)  ::  ind_lu_ev_needle_trees = 4              ! --> ilu_coniferous_forest
+      INTEGER(iwp)  ::  ind_lu_de_needle_trees = 5              ! --> ilu_coniferous_forest
+      INTEGER(iwp)  ::  ind_lu_ev_broad_trees = 6               ! --> ilu_tropical_forest
+      INTEGER(iwp)  ::  ind_lu_de_broad_trees = 7               ! --> ilu_deciduous_forest
+      INTEGER(iwp)  ::  ind_lu_t_grass = 8                      ! --> ilu_grass
+      INTEGER(iwp)  ::  ind_lu_desert = 9                       ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lu_tundra = 10                      ! --> ilu_other
+      INTEGER(iwp)  ::  ind_lu_irr_crops = 11                   ! --> ilu_arable
+      INTEGER(iwp)  ::  ind_lu_semidesert = 12                  ! --> ilu_other
+      INTEGER(iwp)  ::  ind_lu_ice = 13                         ! --> ilu_ice
+      INTEGER(iwp)  ::  ind_lu_marsh = 14                       ! --> ilu_other
+      INTEGER(iwp)  ::  ind_lu_ev_shrubs = 15                   ! --> ilu_mediterrean_scrub
+      INTEGER(iwp)  ::  ind_lu_de_shrubs = 16                   ! --> ilu_mediterrean_scrub
+      INTEGER(iwp)  ::  ind_lu_mixed_forest = 17                ! --> ilu_coniferous_forest (ave(decid+conif))
+      INTEGER(iwp)  ::  ind_lu_intrup_forest = 18               ! --> ilu_other (ave(other+decid))
+
+      ! Water
+      INTEGER(iwp)  ::  ind_luw_user = 0                        ! --> ERROR
+      INTEGER(iwp)  ::  ind_luw_lake = 1                        ! --> ilu_water_inland
+      INTEGER(iwp)  ::  ind_luw_river = 2                       ! --> ilu_water_inland
+      INTEGER(iwp)  ::  ind_luw_ocean = 3                       ! --> ilu_water_sea
+      INTEGER(iwp)  ::  ind_luw_pond = 4                        ! --> ilu_water_inland
+      INTEGER(iwp)  ::  ind_luw_fountain = 5                    ! --> ilu_water_inland 
+
+      ! Pavement
+      INTEGER(iwp)  ::  ind_lup_user = 0                        ! --> ERROR
+      INTEGER(iwp)  ::  ind_lup_asph_conc = 1                   ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_asph = 2                        ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_conc = 3                        ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_sett = 4                        ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_pav_stones = 5                  ! --> ilu_desert 
+      INTEGER(iwp)  ::  ind_lup_cobblest = 6                    ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_metal = 7                       ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_wood = 8                        ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_gravel = 9                      ! --> ilu_desert 
+      INTEGER(iwp)  ::  ind_lup_f_gravel = 10                   ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_pebblest = 11                   ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_woodchips = 12                  ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_tartan = 13                     ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_art_turf = 14                   ! --> ilu_desert
+      INTEGER(iwp)  ::  ind_lup_clay = 15                       ! --> ilu_desert
+
+
+
+      !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
+      INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
+      INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
+      INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
+
+      INTEGER(iwp) ::  part_type
+
+      INTEGER(iwp) :: nwet                ! wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+
+      REAL(kind=wp)  :: dt_chem
+      REAL(kind=wp)  :: dh                
+      REAL(kind=wp)  :: inv_dh
+      REAL(kind=wp)  :: dt_dh
+
+      REAL(kind=wp)  :: dens              ! Density at lowest layer at i,j  
+      REAL(kind=wp)  :: r_aero_lu         ! aerodynamic resistance (s/m) at current surface element
+      REAL(kind=wp)  :: ustar_lu          ! ustar at current surface element
+      REAL(kind=wp)  :: z0h_lu            ! roughness length for heat at current surface element 
+      REAL(kind=wp)  :: glrad             ! rad_sw_in at current surface element 
+      REAL(kind=wp)  :: ppm_to_ugm3       ! conversion factor
+      REAL(kind=wp)  :: relh              ! relative humidity at current surface element
+      REAL(kind=wp)  :: lai               ! leaf area index at current surface element
+      REAL(kind=wp)  :: sai               ! surface area index at current surface element assumed to be lai + 1
+
+      REAL(kind=wp)  :: slinnfac        
+      REAL(kind=wp)  :: visc              ! Viscosity
+      REAL(kind=wp)  :: vs                ! Sedimentation velocity
+      REAL(kind=wp)  :: vd_lu             ! deposition velocity (m/s)
+      REAL(kind=wp)  :: Rs                ! Sedimentaion resistance (s/m)
+      REAL(kind=wp)  :: Rb                ! quasi-laminar boundary layer resistance (s/m)
+      REAL(kind=wp)  :: Rc_tot            ! total canopy resistance Rc (s/m)
+
+      REAL(kind=wp)  :: catm              ! gasconc. at i,j,k=1 in ug/m3
+      REAL(kind=wp)  :: diffc             ! Diffusivity
+
+
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lu       ! budget for LSM vegetation type at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lup      ! budget for LSM pavement type at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_luw      ! budget for LSM water type at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_luu       ! budget for USM walls/roofs at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lug      ! budget for USM green surfaces at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now_lud      ! budget for USM windows at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: bud_now          ! overall budget at current surface element
+      REAL(kind=wp),DIMENSION(nspec)            :: cc               ! concentration i,j,k
+      REAL(kind=wp),DIMENSION(nspec)            :: ccomp_tot        ! total compensation point (ug/m3) (= 0 for species that don't have a compensation)
+      ! For now kept to zero for all species!
+
+      REAL(kind=wp)                             :: ttemp          ! temperatur at i,j,k
+      REAL(kind=wp)                             :: ts             ! surface temperatur in degrees celsius
+      REAL(kind=wp)                             :: qv_tmp         ! surface mixing ratio at current surface element
+
+
+      !-- Particle parameters (PM10 (1), PM25 (2))
+      !-- partsize (diameter in m) ,    rhopart (density in kg/m3),     slipcor (slip correction factor DIMENSIONless, Seinfeld and Pandis 2006, Table 9.3)
+      REAL(wp), DIMENSION(1:3,1:2), PARAMETER   :: particle_pars = RESHAPE( (/ &
+           8.0e-6_wp, 1.14e3_wp, 1.016_wp, & !  1
+           0.7e-6_wp, 1.14e3_wp, 1.082_wp &  !  2
+           /), (/ 3, 2 /) )
+
+
+      LOGICAL                                   ::  match_lsm     ! flag indicating natural-type surface
+      LOGICAL                                   ::  match_usm     ! flag indicating urban-type surface
+
+
+
+      !-- List of names of possible tracers
+      CHARACTER(len=*), PARAMETER  ::  pspecnames(69) = (/ &
+           'NO2           ', &    ! NO2
+           'NO            ', &    ! NO
+           'O3            ', &    ! O3
+           'CO            ', &    ! CO
+           'form          ', &    ! FORM
+           'ald           ', &    ! ALD
+           'pan           ', &    ! PAN
+           'mgly          ', &    ! MGLY
+           'par           ', &    ! PAR
+           'ole           ', &    ! OLE
+           'eth           ', &    ! ETH
+           'tol           ', &    ! TOL
+           'cres          ', &    ! CRES
+           'xyl           ', &    ! XYL
+           'SO4a_f        ', &    ! SO4a_f
+           'SO2           ', &    ! SO2
+           'HNO2          ', &    ! HNO2
+           'CH4           ', &    ! CH4
+           'NH3           ', &    ! NH3
+           'NO3           ', &    ! NO3
+           'OH            ', &    ! OH
+           'HO2           ', &    ! HO2
+           'N2O5          ', &    ! N2O5
+           'SO4a_c        ', &    ! SO4a_c
+           'NH4a_f        ', &    ! NH4a_f
+           'NO3a_f        ', &    ! NO3a_f
+           'NO3a_c        ', &    ! NO3a_c
+           'C2O3          ', &    ! C2O3
+           'XO2           ', &    ! XO2
+           'XO2N          ', &    ! XO2N
+           'cro           ', &    ! CRO
+           'HNO3          ', &    ! HNO3
+           'H2O2          ', &    ! H2O2
+           'iso           ', &    ! ISO
+           'ispd          ', &    ! ISPD
+           'to2           ', &    ! TO2
+           'open          ', &    ! OPEN
+           'terp          ', &    ! TERP
+           'ec_f          ', &    ! EC_f
+           'ec_c          ', &    ! EC_c
+           'pom_f         ', &    ! POM_f
+           'pom_c         ', &    ! POM_c
+           'ppm_f         ', &    ! PPM_f
+           'ppm_c         ', &    ! PPM_c
+           'na_ff         ', &    ! Na_ff
+           'na_f          ', &    ! Na_f
+           'na_c          ', &    ! Na_c
+           'na_cc         ', &    ! Na_cc
+           'na_ccc        ', &    ! Na_ccc
+           'dust_ff       ', &    ! dust_ff
+           'dust_f        ', &    ! dust_f
+           'dust_c        ', &    ! dust_c
+           'dust_cc       ', &    ! dust_cc
+           'dust_ccc      ', &    ! dust_ccc
+           'tpm10         ', &    ! tpm10
+           'tpm25         ', &    ! tpm25
+           'tss           ', &    ! tss
+           'tdust         ', &    ! tdust
+           'tc            ', &    ! tc
+           'tcg           ', &    ! tcg
+           'tsoa          ', &    ! tsoa
+           'tnmvoc        ', &    ! tnmvoc
+           'SOxa          ', &    ! SOxa
+           'NOya          ', &    ! NOya
+           'NHxa          ', &    ! NHxa
+           'NO2_obs       ', &    ! NO2_obs
+           'tpm10_biascorr', &    ! tpm10_biascorr
+           'tpm25_biascorr', &    ! tpm25_biascorr
+           'O3_biascorr   ' /)    ! o3_biascorr
+
+
+
+      ! tracer mole mass:
+      REAL, PARAMETER  ::  specmolm(69) = (/ &
+           xm_O * 2 + xm_N, &                         ! NO2
+           xm_O + xm_N, &                             ! NO
+           xm_O * 3, &                                ! O3
+           xm_C + xm_O, &                             ! CO
+           xm_H * 2 + xm_C + xm_O, &                  ! FORM
+           xm_H * 3 + xm_C * 2 + xm_O, &              ! ALD
+           xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   ! PAN
+           xm_H * 4 + xm_C * 3 + xm_O * 2, &          ! MGLY
+           xm_H * 3 + xm_C, &                         ! PAR
+           xm_H * 3 + xm_C * 2, &                     ! OLE
+           xm_H * 4 + xm_C * 2, &                     ! ETH
+           xm_H * 8 + xm_C * 7, &                     ! TOL
+           xm_H * 8 + xm_C * 7 + xm_O, &              ! CRES
+           xm_H * 10 + xm_C * 8, &                    ! XYL
+           xm_S + xm_O * 4, &                         ! SO4a_f
+           xm_S + xm_O * 2, &                         ! SO2
+           xm_H + xm_O * 2 + xm_N, &                  ! HNO2
+           xm_H * 4 + xm_C, &                         ! CH4
+           xm_H * 3 + xm_N, &                         ! NH3
+           xm_O * 3 + xm_N, &                         ! NO3
+           xm_H + xm_O, &                             ! OH
+           xm_H + xm_O * 2, &                         ! HO2
+           xm_O * 5 + xm_N * 2, &                     ! N2O5
+           xm_S + xm_O * 4, &                         ! SO4a_c
+           xm_H * 4 + xm_N, &                         ! NH4a_f
+           xm_O * 3 + xm_N, &                         ! NO3a_f
+           xm_O * 3 + xm_N, &                         ! NO3a_c
+           xm_C * 2 + xm_O * 3, &                     ! C2O3
+           xm_dummy, &                                ! XO2
+           xm_dummy, &                                ! XO2N
+           xm_dummy, &                                ! CRO
+           xm_H + xm_O * 3 + xm_N, &                  ! HNO3
+           xm_H * 2 + xm_O * 2, &                     ! H2O2
+           xm_H * 8 + xm_C * 5, &                     ! ISO
+           xm_dummy, &                                ! ISPD
+           xm_dummy, &                                ! TO2
+           xm_dummy, &                                ! OPEN
+           xm_H * 16 + xm_C * 10, &                   ! TERP
+           xm_dummy, &                                ! EC_f
+           xm_dummy, &                                ! EC_c
+           xm_dummy, &                                ! POM_f
+           xm_dummy, &                                ! POM_c
+           xm_dummy, &                                ! PPM_f
+           xm_dummy, &                                ! PPM_c
+           xm_Na, &                                   ! Na_ff
+           xm_Na, &                                   ! Na_f
+           xm_Na, &                                   ! Na_c
+           xm_Na, &                                   ! Na_cc
+           xm_Na, &                                   ! Na_ccc
+           xm_dummy, &                                ! dust_ff
+           xm_dummy, &                                ! dust_f
+           xm_dummy, &                                ! dust_c
+           xm_dummy, &                                ! dust_cc
+           xm_dummy, &                                ! dust_ccc
+           xm_dummy, &                                ! tpm10
+           xm_dummy, &                                ! tpm25
+           xm_dummy, &                                ! tss
+           xm_dummy, &                                ! tdust
+           xm_dummy, &                                ! tc
+           xm_dummy, &                                ! tcg
+           xm_dummy, &                                ! tsoa
+           xm_dummy, &                                ! tnmvoc
+           xm_dummy, &                                ! SOxa
+           xm_dummy, &                                ! NOya
+           xm_dummy, &                                ! NHxa
+           xm_O * 2 + xm_N, &                         ! NO2_obs
+           xm_dummy, &                                ! tpm10_biascorr
+           xm_dummy, &                                ! tpm25_biascorr
+           xm_O * 3 /)                                ! o3_biascorr
+
+
+
+      ! Initialize surface element m
+      m = 0
+      k = 0
+      ! LSM or USM surface present at i,j:
+      ! Default surfaces are NOT considered for deposition
+      match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
+      match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
+
+
+      !!!!!!!
+      ! LSM !
+      !!!!!!!
+
+      IF ( match_lsm )  THEN
+
+         ! Get surface element information at i,j:
+         m = surf_lsm_h%start_index(j,i)
+
+         ! Get corresponding k
+         k = surf_lsm_h%k(m)
+
+         ! Get needed variables for surface element m
+         ustar_lu  = surf_lsm_h%us(m)
+         z0h_lu    = surf_lsm_h%z0h(m)
+         r_aero_lu = surf_lsm_h%r_a(m)
+         glrad     = surf_lsm_h%rad_sw_in(m)
+         lai = surf_lsm_h%lai(m)
+         sai = lai + 1
+
+         ! For small grid spacing neglect R_a
+         IF (dzw(k) <= 1.0) THEN
+            r_aero_lu = 0.0_wp
+         ENDIF
+
+         !Initialize lu's
+         lu_palm = 0
+         lu_dep = 0
+         lup_palm = 0
+         lup_dep = 0
+         luw_palm = 0
+         luw_dep = 0
+
+         !Initialize budgets
+         bud_now_lu  = 0.0_wp
+         bud_now_lup = 0.0_wp
+         bud_now_luw = 0.0_wp
+
+
+         ! Get land use for i,j and assign to DEPAC lu
+         IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN
+            lu_palm = surf_lsm_h%vegetation_type(m)
+            IF (lu_palm == ind_lu_user) THEN
+               message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
+               CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
+            ELSEIF (lu_palm == ind_lu_b_soil) THEN
+               lu_dep = 9
+            ELSEIF (lu_palm == ind_lu_mixed_crops) THEN
+               lu_dep = 2
+            ELSEIF (lu_palm == ind_lu_s_grass) THEN
+               lu_dep = 1
+            ELSEIF (lu_palm == ind_lu_ev_needle_trees) THEN
+               lu_dep = 4
+            ELSEIF (lu_palm == ind_lu_de_needle_trees) THEN
+               lu_dep = 4
+            ELSEIF (lu_palm == ind_lu_ev_broad_trees) THEN
+               lu_dep = 12
+            ELSEIF (lu_palm == ind_lu_de_broad_trees) THEN
+               lu_dep = 5
+            ELSEIF (lu_palm == ind_lu_t_grass) THEN
+               lu_dep = 1
+            ELSEIF (lu_palm == ind_lu_desert) THEN
+               lu_dep = 9
+            ELSEIF (lu_palm == ind_lu_tundra ) THEN
+               lu_dep = 8
+            ELSEIF (lu_palm == ind_lu_irr_crops) THEN
+               lu_dep = 2
+            ELSEIF (lu_palm == ind_lu_semidesert) THEN
+               lu_dep = 8
+            ELSEIF (lu_palm == ind_lu_ice) THEN
+               lu_dep = 10
+            ELSEIF (lu_palm == ind_lu_marsh) THEN
+               lu_dep = 8
+            ELSEIF (lu_palm == ind_lu_ev_shrubs) THEN
+               lu_dep = 14
+            ELSEIF (lu_palm == ind_lu_de_shrubs ) THEN
+               lu_dep = 14
+            ELSEIF (lu_palm == ind_lu_mixed_forest) THEN
+               lu_dep = 4
+            ELSEIF (lu_palm == ind_lu_intrup_forest) THEN
+               lu_dep = 8      
+            ENDIF
+         ENDIF
+
+         IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN
+            lup_palm = surf_lsm_h%pavement_type(m)
+            IF (lup_palm == ind_lup_user) THEN
+               message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
+               CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
+            ELSEIF (lup_palm == ind_lup_asph_conc) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_asph) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm ==  ind_lup_conc) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm ==  ind_lup_sett) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_pav_stones) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_cobblest) THEN
+               lup_dep = 9       
+            ELSEIF (lup_palm == ind_lup_metal) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_wood) THEN
+               lup_dep = 9    
+            ELSEIF (lup_palm == ind_lup_gravel) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_f_gravel) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_pebblest) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_woodchips) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_tartan) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_art_turf) THEN
+               lup_dep = 9
+            ELSEIF (lup_palm == ind_lup_clay) THEN
+               lup_dep = 9
+            ENDIF
+         ENDIF
+
+         IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN
+            luw_palm = surf_lsm_h%water_type(m)     
+            IF (luw_palm == ind_luw_user) THEN
+               message_string = 'No water type defined. Please define water type to enable deposition calculation'
+               CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
+            ELSEIF (luw_palm ==  ind_luw_lake) THEN
+               luw_dep = 13
+            ELSEIF (luw_palm == ind_luw_river) THEN
+               luw_dep = 13
+            ELSEIF (luw_palm == ind_luw_ocean) THEN
+               luw_dep = 6
+            ELSEIF (luw_palm == ind_luw_pond) THEN
+               luw_dep = 13
+            ELSEIF (luw_palm == ind_luw_fountain) THEN
+               luw_dep = 13 
+            ENDIF
+         ENDIF
+
+
+
+         ! Set wetness indicator to dry or wet for lsm vegetation or pavement
+         IF (surf_lsm_h%c_liq(m) > 0) THEN
+            nwet = 1
+         ELSE
+            nwet = 0
+         ENDIF
+
+         ! Compute length of time step 
+         IF ( call_chem_at_all_substeps )  THEN
+            dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
+         ELSE
+            dt_chem = dt_3d
+         ENDIF
+
+
+         dh = dzw(k)
+         inv_dh = 1.0_wp / dh
+         dt_dh = dt_chem/dh
+
+         ! Concentration at i,j,k
+         DO lsp = 1, nspec
+            cc(lsp) = chem_species(lsp)%conc(k,j,i)
+         ENDDO
+
+
+         ! Temperature column at i,j
+         ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
+
+         ! Surface temperature in degrees Celcius:
+         ts       = ttemp - 273.15 ! in degrees celcius
+
+         ! Viscosity of air in lowest layer
+         visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
+
+         ! Air density in lowest layer
+         dens = rho_air_zw(k)
+
+         ! Calculate relative humidity from specific humidity for DEPAC
+         qv_tmp = max(q(k,j,i),0.0_wp)
+         relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )
+
+
+
+         ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
+         ! for each surface fraction. Then derive overall budget taking into account the surface fractions.
+
+         ! Vegetation
+         IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN
+
+
+            slinnfac = 1.0_wp
+
+            ! Get vd
+            DO lsp = 1, nvar
+               !Initialize 
+               vs = 0.0_wp
+               vd_lu = 0.0_wp
+               Rs = 0.0_wp
+               Rb = 0.0_wp
+               Rc_tot = 0.0_wp
+               IF (spc_names(lsp) == 'PM10') THEN
+                  part_type = 1
+                  ! sedimentation velocity in lowest layer
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lu_dep,  &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lu(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSEIF (spc_names(lsp) == 'PM25') THEN
+                  part_type = 2
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lu_dep , &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lu(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSE
+
+                  ! GASES
+
+                  ! Read spc_name of current species for gas parameter
+
+                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
+                     i_pspec = 0
+                     DO pspec = 1, 69
+                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
+                           i_pspec = pspec
+                        END IF
+                     ENDDO
+
+                  ELSE
+                     ! Species for now not deposited
+                     CYCLE
+                  ENDIF
+
+                  ! Factor used for conversion from ppb to ug/m3 :
+                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+                  !    c           1e-9              xm_tracer         1e9       /       xm_air            dens
+                  ! thus:
+                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+                  ! Use density at lowest layer:
+
+                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
+
+                  ! atmospheric concentration in DEPAC is requested in ug/m3:
+                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+
+                  ! Diffusivity for DEPAC relevant gases
+                  ! Use default value 
+                  diffc            = 0.11e-4
+                  ! overwrite with known coefficients of diffusivity from Massman (1998)
+                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
+                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
+                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
+                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
+                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
+                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
+                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
+
+
+                  ! Get quasi-laminar boundary layer resistance Rb:
+                  CALL get_rb_cell( (lu_dep == ilu_water_sea) .or. (lu_dep == ilu_water_inland), &
+                       z0h_lu, ustar_lu, diffc, &
+                       Rb)
+
+                  ! Get Rc_tot
+                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                       relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
+                       r_aero_lu , Rb)
+
+
+                  ! Calculate budget
+                  IF (Rc_tot .le. 0.0) THEN
+
+                     bud_now_lu(lsp) = 0.0_wp
+
+                  ELSE
+
+                     ! Compute exchange velocity for current lu:
+                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
+
+                     bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                          (1.0 - exp(-vd_lu*dt_dh ))*dh
+                  ENDIF
+
+               ENDIF
+            ENDDO
+         ENDIF
+
+
+         ! Pavement
+         IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN
+
+
+            ! No vegetation on pavements:
+            lai = 0.0_wp
+            sai = 0.0_wp
+            
+            slinnfac = 1.0_wp
+
+            ! Get vd
+            DO lsp = 1, nvar
+               !Initialize
+               vs = 0.0_wp
+               vd_lu = 0.0_wp
+               Rs = 0.0_wp
+               Rb = 0.0_wp
+               Rc_tot = 0.0_wp
+               IF (spc_names(lsp) == 'PM10') THEN
+                  part_type = 1
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lup_dep,  &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lup(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSEIF (spc_names(lsp) == 'PM25') THEN
+                  part_type = 2
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lup_dep, &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lup(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSE
+
+                  ! GASES
+
+                  ! Read spc_name of current species for gas parameter
+
+                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
+                     i_pspec = 0
+                     DO pspec = 1, 69
+                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
+                           i_pspec = pspec
+                        END IF
+                     ENDDO
+
+                  ELSE
+                     ! Species for now not deposited
+                     CYCLE
+                  ENDIF
+
+                  ! Factor used for conversion from ppb to ug/m3 :
+                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
+                  ! thus:
+                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+                  ! Use density at lowest layer:
+
+                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
+
+                  ! atmospheric concentration in DEPAC is requested in ug/m3:
+                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+
+                  ! Diffusivity for DEPAC relevant gases
+                  ! Use default value 
+                  diffc            = 0.11e-4
+                  ! overwrite with known coefficients of diffusivity from Massman (1998)
+                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
+                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
+                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
+                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
+                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
+                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
+                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
+
+
+                  ! Get quasi-laminar boundary layer resistance Rb:
+                  CALL get_rb_cell( (lup_dep == ilu_water_sea) .or. (lup_dep == ilu_water_inland), &
+                       z0h_lu, ustar_lu, diffc, &
+                       Rb)
+
+
+                  ! Get Rc_tot
+                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                       relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
+                       r_aero_lu , Rb)
+
+
+                  ! Calculate budget
+                  IF (Rc_tot .le. 0.0) THEN
+
+                     bud_now_lup(lsp) = 0.0_wp
+
+                  ELSE
+
+                     ! Compute exchange velocity for current lu:
+                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
+
+                     bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                          (1.0 - exp(-vd_lu*dt_dh ))*dh
+                  ENDIF
+
+
+               ENDIF
+            ENDDO
+         ENDIF
+
+
+         ! Water
+         IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN
+
+
+            ! No vegetation on water:
+            lai = 0.0_wp
+            sai = 0.0_wp
+            
+            slinnfac = 1.0_wp
+
+            ! Get vd
+            DO lsp = 1, nvar
+               !Initialize
+               vs = 0.0_wp
+               vd_lu = 0.0_wp
+               Rs = 0.0_wp
+               Rb = 0.0_wp
+               Rc_tot = 0.0_wp 
+               IF (spc_names(lsp) == 'PM10') THEN
+                  part_type = 1
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       luw_dep, &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_luw(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSEIF (spc_names(lsp) == 'PM25') THEN
+                  part_type = 2
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       luw_dep, &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_luw(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSE
+
+                  ! GASES
+
+                  ! Read spc_name of current species for gas PARAMETER
+
+                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
+                     i_pspec = 0
+                     DO pspec = 1, 69
+                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
+                           i_pspec = pspec
+                        END IF
+                     ENDDO
+
+                  ELSE
+                     ! Species for now not deposited
+                     CYCLE
+                  ENDIF
+
+                  ! Factor used for conversion from ppb to ug/m3 :
+                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
+                  ! thus:
+                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+                  ! Use density at lowest layer:
+
+                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
+
+                  ! atmospheric concentration in DEPAC is requested in ug/m3:
+                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+
+                  ! Diffusivity for DEPAC relevant gases
+                  ! Use default value 
+                  diffc            = 0.11e-4
+                  ! overwrite with known coefficients of diffusivity from Massman (1998)
+                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
+                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
+                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
+                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
+                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
+                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
+                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
+
+
+                  ! Get quasi-laminar boundary layer resistance Rb:
+                  CALL get_rb_cell( (luw_dep == ilu_water_sea) .or. (luw_dep == ilu_water_inland), &
+                       z0h_lu, ustar_lu, diffc, &
+                       Rb)
+
+                  ! Get Rc_tot
+                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                       relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
+                       r_aero_lu , Rb)
+
+
+                  ! Calculate budget 
+                  IF (Rc_tot .le. 0.0) THEN
+
+                     bud_now_luw(lsp) = 0.0_wp
+
+                  ELSE
+
+                     ! Compute exchange velocity for current lu:
+                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
+
+                     bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                          (1.0 - exp(-vd_lu*dt_dh ))*dh
+                  ENDIF
+
+               ENDIF
+            ENDDO
+         ENDIF
+
+
+         bud_now = 0.0_wp
+         ! Calculate budget for surface m and adapt concentration
+         DO lsp = 1, nspec
+
+
+            bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m)*bud_now_lu(lsp) + &
+                 surf_lsm_h%frac(ind_pav_green,m)*bud_now_lup(lsp) + &
+                 surf_lsm_h%frac(ind_wat_win,m)*bud_now_luw(lsp)
+
+            ! Compute new concentration, add contribution of all exchange processes:
+            cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh
+            
+            chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp))
+ 
+         ENDDO
+
+      ENDIF
+
+
+
+      !!!!!!!
+      ! USM !
+      !!!!!!!    
+
+      IF ( match_usm )  THEN
+
+         ! Get surface element information at i,j:
+         m = surf_usm_h%start_index(j,i)
+
+         k = surf_usm_h%k(m)
+
+         ! Get needed variables for surface element m
+         ustar_lu  = surf_usm_h%us(m)
+         z0h_lu    = surf_usm_h%z0h(m)
+         r_aero_lu = surf_usm_h%r_a(m)
+         glrad     = surf_usm_h%rad_sw_in(m)
+         lai = surf_usm_h%lai(m)
+         sai = lai + 1
+
+         ! For small grid spacing neglect R_a
+         IF (dzw(k) <= 1.0) THEN
+            r_aero_lu = 0.0_wp
+         ENDIF
+
+         !Initialize lu's
+         luu_palm = 0
+         luu_dep = 0
+         lug_palm = 0
+         lug_dep = 0
+         lud_palm = 0
+         lud_dep = 0
+
+         !Initialize budgets
+         bud_now_luu  = 0.0_wp
+         bud_now_lug = 0.0_wp
+         bud_now_lud = 0.0_wp
+
+
+         ! Get land use for i,j and assign to DEPAC lu
+         IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN
+            ! For green urban surfaces (e.g. green roofs
+            ! assume LU short grass
+            lug_palm = ind_lu_s_grass
+            IF (lug_palm == ind_lu_user) THEN
+               message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
+               CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
+            ELSEIF (lug_palm == ind_lu_b_soil) THEN
+               lug_dep = 9
+            ELSEIF (lug_palm == ind_lu_mixed_crops) THEN
+               lug_dep = 2
+            ELSEIF (lug_palm == ind_lu_s_grass) THEN
+               lug_dep = 1
+            ELSEIF (lug_palm == ind_lu_ev_needle_trees) THEN
+               lug_dep = 4
+            ELSEIF (lug_palm == ind_lu_de_needle_trees) THEN
+               lug_dep = 4
+            ELSEIF (lug_palm == ind_lu_ev_broad_trees) THEN
+               lug_dep = 12
+            ELSEIF (lug_palm == ind_lu_de_broad_trees) THEN
+               lug_dep = 5
+            ELSEIF (lug_palm == ind_lu_t_grass) THEN
+               lug_dep = 1
+            ELSEIF (lug_palm == ind_lu_desert) THEN
+               lug_dep = 9
+            ELSEIF (lug_palm == ind_lu_tundra ) THEN
+               lug_dep = 8
+            ELSEIF (lug_palm == ind_lu_irr_crops) THEN
+               lug_dep = 2
+            ELSEIF (lug_palm == ind_lu_semidesert) THEN
+               lug_dep = 8
+            ELSEIF (lug_palm == ind_lu_ice) THEN
+               lug_dep = 10
+            ELSEIF (lug_palm == ind_lu_marsh) THEN
+               lug_dep = 8
+            ELSEIF (lug_palm == ind_lu_ev_shrubs) THEN
+               lug_dep = 14
+            ELSEIF (lug_palm == ind_lu_de_shrubs ) THEN
+               lug_dep = 14
+            ELSEIF (lug_palm == ind_lu_mixed_forest) THEN
+               lug_dep = 4
+            ELSEIF (lug_palm == ind_lu_intrup_forest) THEN
+               lug_dep = 8      
+            ENDIF
+         ENDIF
+
+         IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN
+            ! For walls in USM assume concrete walls/roofs,
+            ! assumed LU class desert as also assumed for
+            ! pavements in LSM
+            luu_palm = ind_lup_conc
+            IF (luu_palm == ind_lup_user) THEN
+               message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
+               CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
+            ELSEIF (luu_palm == ind_lup_asph_conc) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_asph) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm ==  ind_lup_conc) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm ==  ind_lup_sett) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_pav_stones) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_cobblest) THEN
+               luu_dep = 9       
+            ELSEIF (luu_palm == ind_lup_metal) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_wood) THEN
+               luu_dep = 9    
+            ELSEIF (luu_palm == ind_lup_gravel) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_f_gravel) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_pebblest) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_woodchips) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_tartan) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_art_turf) THEN
+               luu_dep = 9
+            ELSEIF (luu_palm == ind_lup_clay) THEN
+               luu_dep = 9
+            ENDIF
+         ENDIF
+
+         IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN
+            ! For windows in USM assume metal as this is
+            ! as close as we get, assumed LU class desert
+            ! as also assumed for pavements in LSM
+            lud_palm = ind_lup_metal     
+            IF (lud_palm == ind_lup_user) THEN
+               message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
+               CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
+            ELSEIF (lud_palm == ind_lup_asph_conc) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_asph) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm ==  ind_lup_conc) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm ==  ind_lup_sett) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_pav_stones) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_cobblest) THEN
+               lud_dep = 9       
+            ELSEIF (lud_palm == ind_lup_metal) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_wood) THEN
+               lud_dep = 9    
+            ELSEIF (lud_palm == ind_lup_gravel) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_f_gravel) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_pebblest) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_woodchips) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_tartan) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_art_turf) THEN
+               lud_dep = 9
+            ELSEIF (lud_palm == ind_lup_clay) THEN
+               lud_dep = 9
+            ENDIF
+         ENDIF
+
+
+         ! TODO: Activate these lines as soon as new ebsolver branch is merged:
+         ! Set wetness indicator to dry or wet for usm vegetation or pavement
+         !IF (surf_usm_h%c_liq(m) > 0) THEN
+         !   nwet = 1
+         !ELSE
+         nwet = 0
+         !ENDIF
+
+         ! Compute length of time step 
+         IF ( call_chem_at_all_substeps )  THEN
+            dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
+         ELSE
+            dt_chem = dt_3d
+         ENDIF
+
+
+         dh = dzw(k)
+         inv_dh = 1.0_wp / dh
+         dt_dh = dt_chem/dh
+
+         ! Concentration at i,j,k
+         DO lsp = 1, nspec
+            cc(lsp) = chem_species(lsp)%conc(k,j,i)
+         ENDDO
+
+         ! Temperature at i,j,k
+         ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
+
+         ! Surface temperature in degrees Celcius:
+         ts       = ttemp - 273.15 ! in degrees celcius
+
+         ! Viscosity of air in lowest layer
+         visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
+
+         ! Air density in lowest layer
+         dens = rho_air_zw(k)
+
+         ! Calculate relative humidity from specific humidity for DEPAC
+         qv_tmp = max(q(k,j,i),0.0_wp)
+         relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )
+
+
+
+         ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
+         ! for each surface fraction. Then derive overall budget taking into account the surface fractions.
+
+         ! Walls/roofs
+         IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN
+
+
+            ! No vegetation on non-green walls:
+            lai = 0.0_wp
+            sai = 0.0_wp
+            
+            slinnfac = 1.0_wp
+
+            ! Get vd
+            DO lsp = 1, nvar
+               !Initialize 
+               vs = 0.0_wp
+               vd_lu = 0.0_wp
+               Rs = 0.0_wp
+               Rb = 0.0_wp
+               Rc_tot = 0.0_wp
+               IF (spc_names(lsp) == 'PM10') THEN
+                  part_type = 1
+                  ! sedimentation velocity in lowest layer
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       luu_dep,  &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_luu(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSEIF (spc_names(lsp) == 'PM25') THEN
+                  part_type = 2
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       luu_dep , &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_luu(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSE
+
+                  ! GASES
+
+                  ! Read spc_name of current species for gas parameter
+
+                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
+                     i_pspec = 0
+                     DO pspec = 1, 69
+                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
+                           i_pspec = pspec
+                        END IF
+                     ENDDO
+
+                  ELSE
+                     ! Species for now not deposited
+                     CYCLE
+                  ENDIF
+
+                  ! Factor used for conversion from ppb to ug/m3 :
+                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+                  !    c           1e-9              xm_tracer         1e9       /       xm_air            dens
+                  ! thus:
+                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+                  ! Use density at lowest layer:
+
+                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
+
+                  ! atmospheric concentration in DEPAC is requested in ug/m3:
+                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+
+                  ! Diffusivity for DEPAC relevant gases
+                  ! Use default value 
+                  diffc            = 0.11e-4
+                  ! overwrite with known coefficients of diffusivity from Massman (1998)
+                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
+                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
+                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
+                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
+                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
+                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
+                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
+
+
+                  ! Get quasi-laminar boundary layer resistance Rb:
+                  CALL get_rb_cell( (luu_dep == ilu_water_sea) .or. (luu_dep == ilu_water_inland), &
+                       z0h_lu, ustar_lu, diffc, &
+                       Rb)
+
+                  ! Get Rc_tot
+                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                       relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
+                       r_aero_lu , Rb)
+
+
+                  ! Calculate budget
+                  IF (Rc_tot .le. 0.0) THEN
+
+                     bud_now_luu(lsp) = 0.0_wp
+
+                  ELSE
+
+                     ! Compute exchange velocity for current lu:
+                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
+
+                     bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                          (1.0 - exp(-vd_lu*dt_dh ))*dh
+                  ENDIF
+
+               ENDIF
+            ENDDO
+         ENDIF
+
+
+         ! Green usm surfaces
+         IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN
+
+
+            slinnfac = 1.0_wp
+
+            ! Get vd
+            DO lsp = 1, nvar
+               !Initialize
+               vs = 0.0_wp
+               vd_lu = 0.0_wp
+               Rs = 0.0_wp
+               Rb = 0.0_wp
+               Rc_tot = 0.0_wp
+               IF (spc_names(lsp) == 'PM10') THEN
+                  part_type = 1
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lug_dep,  &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lug(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSEIF (spc_names(lsp) == 'PM25') THEN
+                  part_type = 2
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lug_dep, &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lug(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSE
+
+                  ! GASES
+
+                  ! Read spc_name of current species for gas parameter
+
+                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
+                     i_pspec = 0
+                     DO pspec = 1, 69
+                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
+                           i_pspec = pspec
+                        END IF
+                     ENDDO
+
+                  ELSE
+                     ! Species for now not deposited
+                     CYCLE
+                  ENDIF
+
+                  ! Factor used for conversion from ppb to ug/m3 :
+                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
+                  ! thus:
+                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+                  ! Use density at lowest layer:
+
+                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
+
+                  ! atmospheric concentration in DEPAC is requested in ug/m3:
+                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+
+                  ! Diffusivity for DEPAC relevant gases
+                  ! Use default value 
+                  diffc            = 0.11e-4
+                  ! overwrite with known coefficients of diffusivity from Massman (1998)
+                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
+                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
+                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
+                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
+                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
+                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
+                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
+
+
+                  ! Get quasi-laminar boundary layer resistance Rb:
+                  CALL get_rb_cell( (lug_dep == ilu_water_sea) .or. (lug_dep == ilu_water_inland), &
+                       z0h_lu, ustar_lu, diffc, &
+                       Rb)
+
+
+                  ! Get Rc_tot
+                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                       relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
+                       r_aero_lu , Rb)
+
+
+                  ! Calculate budget
+                  IF (Rc_tot .le. 0.0) THEN
+
+                     bud_now_lug(lsp) = 0.0_wp
+
+                  ELSE
+
+                     ! Compute exchange velocity for current lu:
+                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
+
+                     bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                          (1.0 - exp(-vd_lu*dt_dh ))*dh
+                  ENDIF
+
+
+               ENDIF
+            ENDDO
+         ENDIF
+
+
+         ! Windows
+         IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN
+
+
+            ! No vegetation on windows:
+            lai = 0.0_wp
+            sai = 0.0_wp
+            
+            slinnfac = 1.0_wp
+
+            ! Get vd
+            DO lsp = 1, nvar
+               !Initialize
+               vs = 0.0_wp
+               vd_lu = 0.0_wp
+               Rs = 0.0_wp
+               Rb = 0.0_wp
+               Rc_tot = 0.0_wp 
+               IF (spc_names(lsp) == 'PM10') THEN
+                  part_type = 1
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lud_dep, &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lud(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSEIF (spc_names(lsp) == 'PM25') THEN
+                  part_type = 2
+                  ! sedimentation velocity in lowest layer:
+                  vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       visc)
+
+
+                  CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                       vs, &
+                       particle_pars(ind_p_size, part_type), &
+                       particle_pars(ind_p_slip, part_type), &
+                       nwet, ttemp, dens, visc, &
+                       lud_dep, &
+                       r_aero_lu, ustar_lu)
+
+                  bud_now_lud(lsp) = - cc(lsp) * &
+                       (1.0 - exp(-vd_lu*dt_dh ))*dh
+
+
+               ELSE
+
+                  ! GASES
+
+                  ! Read spc_name of current species for gas PARAMETER
+
+                  IF (ANY(pspecnames(:) == spc_names(lsp))) THEN
+                     i_pspec = 0
+                     DO pspec = 1, 69
+                        IF (pspecnames(pspec) == spc_names(lsp)) THEN
+                           i_pspec = pspec
+                        END IF
+                     ENDDO
+
+                  ELSE
+                     ! Species for now not deposited
+                     CYCLE
+                  ENDIF
+
+                  ! Factor used for conversion from ppb to ug/m3 :
+                  !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
+                  !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
+                  !    c           1e-9               xm_tracer         1e9       /       xm_air            dens
+                  ! thus:
+                  !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
+                  ! Use density at lowest layer:
+
+                  ppm_to_ugm3 =  (dens/xm_air)*0.001_wp  ! (mole air)/m3
+
+                  ! atmospheric concentration in DEPAC is requested in ug/m3:
+                  !  ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
+                  catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
+
+                  ! Diffusivity for DEPAC relevant gases
+                  ! Use default value 
+                  diffc            = 0.11e-4
+                  ! overwrite with known coefficients of diffusivity from Massman (1998)
+                  IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
+                  IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
+                  IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
+                  IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
+                  IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
+                  IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
+                  IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
+
+
+                  ! Get quasi-laminar boundary layer resistance Rb:
+                  CALL get_rb_cell( (lud_dep == ilu_water_sea) .or. (lud_dep == ilu_water_inland), &
+                       z0h_lu, ustar_lu, diffc, &
+                       Rb)
+
+                  ! Get Rc_tot
+                  CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                       relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
+                       r_aero_lu , Rb)
+
+
+                  ! Calculate budget 
+                  IF (Rc_tot .le. 0.0) THEN
+
+                     bud_now_lud(lsp) = 0.0_wp
+
+                  ELSE
+
+                     ! Compute exchange velocity for current lu:
+                     vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot )
+
+                     bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                          (1.0 - exp(-vd_lu*dt_dh ))*dh
+                  ENDIF
+
+               ENDIF
+            ENDDO
+         ENDIF
+
+
+         bud_now = 0.0_wp
+         ! Calculate budget for surface m and adapt concentration
+         DO lsp = 1, nspec
+
+
+            bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m)*bud_now_luu(lsp) + &
+                 surf_usm_h%frac(ind_pav_green,m)*bud_now_lug(lsp) + &
+                 surf_usm_h%frac(ind_wat_win,m)*bud_now_lud(lsp)
+
+            ! Compute new concentration, add contribution of all exchange processes:
+            cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh
+
+            chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp))
+
+         ENDDO
+
+      ENDIF
+
+
+
+    END SUBROUTINE chem_depo
+
+  
+
+!----------------------------------------------------------------------------------
+!
+! DEPAC:
+! Code of the DEPAC routine and corresponding subroutines below from the DEPAC
+! module of the LOTOS-EUROS model (Manders et al., 2017)
+!   
+! Original DEPAC routines by RIVM and TNO (2015), for Documentation see
+! van Zanten et al., 2010.
+!---------------------------------------------------------------------------------
+!    
+!----------------------------------------------------------------------------------    
+!
+! depac: compute total canopy (or surface) resistance Rc for gases
+!----------------------------------------------------------------------------------
+
+    SUBROUTINE drydepos_gas_depac(compnam, day_of_year, lat, t, ust, glrad, sinphi, &
+         rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc, &
+         ra, rb)  
+
+
+      ! The last four rows of depac arguments are OPTIONAL:
+      !
+      ! A. compute Rc_tot without compensation points (ccomp_tot will be zero):
+      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
+      !
+      ! B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
+      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
+      !                  c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
+      !
+      ! C. compute effective Rc based on compensation points (used for OPS):
+      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
+      !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
+      !                 ra, rb, rc_eff)
+      ! X1. Extra (OPTIONAL) output variables: 
+      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
+      !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
+      !                 ra, rb, rc_eff, &
+      !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
+      ! X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): 
+      !     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
+      !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
+      !                 ra, rb, rc_eff, &
+      !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
+      !                 calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
+
+
+      IMPLICIT NONE
+
+      CHARACTER(len=*)          , INTENT(in)       :: compnam          ! component name
+      ! 'HNO3','NO','NO2','O3','SO2','NH3'
+      INTEGER(iwp)              , INTENT(in)       :: day_of_year      ! day of year, 1 ... 365 (366)
+      REAL(kind=wp)             , INTENT(in)       :: lat                   ! latitude Northern hemisphere (degrees) (DEPAC cannot be used for S. hemisphere)
+      REAL(kind=wp)            , INTENT(in)        :: t                ! temperature (C)
+      ! NB discussion issue is temp T_2m or T_surf or T_leaf?
+      REAL(kind=wp)             , INTENT(in)       :: ust              ! friction velocity (m/s)
+      REAL(kind=wp)             , INTENT(in)       :: glrad            ! global radiation (W/m2)
+      REAL(kind=wp)             , INTENT(in)       :: sinphi           ! sin of solar elevation angle
+      REAL(kind=wp)             , INTENT(in)       :: rh               ! relative humidity (%)
+      REAL(kind=wp)             , INTENT(in)       :: lai              ! one-sidedleaf area index (-)
+      REAL(kind=wp)             , INTENT(in)       :: sai              ! surface area index (-) (lai + branches and stems)
+      REAL(kind=wp)             , INTENT(in)       :: diffc            ! Diffusivity
+      INTEGER(iwp)              , INTENT(in)           :: nwet             ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+      INTEGER(iwp)              , INTENT(in)       :: lu               ! land use type, lu = 1,...,nlu
+      INTEGER(iwp)              , INTENT(in)       :: iratns           ! index for NH3/SO2 ratio used for SO2;
+      ! iratns = 1: low NH3/SO2
+      ! iratns = 2: high NH3/SO2
+      ! iratns = 3: very low NH3/SO2
+      REAL(kind=wp)             , INTENT(out)      :: rc_tot           ! total canopy resistance Rc (s/m)
+      REAL(kind=wp)             , INTENT(out)      :: ccomp_tot        ! total compensation point (ug/m3) [= 0 for species that don't have a compensation 
+      REAL(kind=wp)             ,INTENT(in)        :: p                ! pressure (Pa)
+      REAL(kind=wp), INTENT(in)  :: catm            ! actual atmospheric concentration (ug/m3)
+      REAL(kind=wp), INTENT(in)  :: ra              ! aerodynamic resistance (s/m)
+      REAL(kind=wp), INTENT(in)  :: rb              ! boundary layer resistance (s/m)
+
+      ! Local variables:
+      REAL(kind=wp)                           :: laimax          ! maximum leaf area index (-)
+      LOGICAL                        :: ready           ! Rc has been set
+      ! = 1 -> constant Rc
+      ! = 2 -> temperature dependent Rc
+      REAL(kind=wp)                           :: gw              ! external leaf conductance (m/s)
+      REAL(kind=wp)                           :: gstom           ! stomatal conductance (m/s)
+      REAL(kind=wp)                           :: gsoil_eff       ! effective soil conductance (m/s)
+      REAL(kind=wp)                           :: gc_tot          ! total canopy conductance (m/s)
+      REAL(kind=wp)                           :: cw              ! external leaf surface compensation point (ug/m3)
+      REAL(kind=wp)                           :: cstom           ! stomatal compensation point (ug/m3)
+      REAL(kind=wp)                           :: csoil           ! soil compensation point (ug/m3)
+
+      ! Vegetation indicators:
+      LOGICAL            :: LAI_present ! leaves are present for current land use type
+      LOGICAL            :: SAI_present ! vegetation is present for current land use type
+
+      ! Component number taken from component name, paramteres matched with include files
+      INTEGER(iwp)           ::  icmp
+
+      ! component numbers:
+      INTEGER(iwp), PARAMETER  ::  icmp_o3   = 1
+      INTEGER(iwp), PARAMETER  ::  icmp_so2  = 2
+      INTEGER(iwp), PARAMETER  ::  icmp_no2  = 3
+      INTEGER(iwp), PARAMETER  ::  icmp_no   = 4
+      INTEGER(iwp), PARAMETER  ::  icmp_nh3  = 5
+      INTEGER(iwp), PARAMETER  ::  icmp_co   = 6
+      INTEGER(iwp), PARAMETER  ::  icmp_no3  = 7
+      INTEGER(iwp), PARAMETER  ::  icmp_hno3 = 8
+      INTEGER(iwp), PARAMETER  ::  icmp_n2o5 = 9
+      INTEGER(iwp), PARAMETER  ::  icmp_h2o2 = 10
+
+
+
+      ! Define component number
+      SELECT CASE ( trim(compnam) )
+
+      CASE ( 'O3', 'o3' )
+         icmp = icmp_o3
+
+      CASE ( 'SO2', 'so2' )
+         icmp = icmp_so2
+
+      CASE ( 'NO2', 'no2' )
+         icmp = icmp_no2
+
+      CASE ( 'NO', 'no' )
+         icmp = icmp_no 
+
+      CASE ( 'NH3', 'nh3' )
+         icmp = icmp_nh3
+
+      CASE ( 'CO', 'co' )
+         icmp = icmp_co
+
+      CASE ( 'NO3', 'no3' )
+         icmp = icmp_no3
+
+      CASE ( 'HNO3', 'hno3' )
+         icmp = icmp_hno3
+
+      CASE ( 'N2O5', 'n2o5' )
+         icmp = icmp_n2o5
+
+      CASE ( 'H2O2', 'h2o2' )
+         icmp = icmp_h2o2
+
+      CASE default
+         !Component not part of DEPAC --> not deposited
+         RETURN
+
+      END SELECT
+
+      ! inititalize
+      gw        = 0.
+      gstom     = 0.
+      gsoil_eff = 0.
+      gc_tot    = 0.
+      cw        = 0.
+      cstom     = 0.
+      csoil     = 0.
+
+
+      ! Check whether vegetation is present (in that CASE the leaf or surface area index > 0):
+      LAI_present = (lai .gt. 0.0)
+      SAI_present = (sai .gt. 0.0)
+
+      ! Set Rc (i.e. rc_tot) in special cases:
+      CALL rc_special(icmp,compnam,lu,t, nwet,rc_tot,ready,ccomp_tot)
+
+      ! set effective resistance equal to total resistance, this will be changed for compensation points
+      !IF ( present(rc_eff) ) then
+      !   rc_eff = rc_tot
+      !END IF
+
+      ! IF Rc is not set:
+      IF (.not. ready) then
+
+         ! External conductance:
+         CALL rc_gw(compnam,iratns,t,rh,nwet,SAI_present,sai,gw)         
+
+         ! Stomatal conductance:
+         ! CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p,&
+         !     smi=smi,calc_stom_o3flux=calc_stom_o3flux)
+         CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p)
+         ! Effective soil conductance:
+         CALL rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff)
+
+         ! Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
+         CALL rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot)
+
+         ! Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
+         ! CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
+         !     LAI_present, SAI_present, &
+         !     ccomp_tot, &
+         !     catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
+         !     c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
+         !     tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
+
+         ! Effective Rc based on compensation points:
+         ! IF ( present(rc_eff) ) then
+         ! check on required arguments:
+         !IF ( (.not. present(catm)) .or. (.not. present(ra)) .or. (.not. present(rb)) ) then
+         !   stop 'output argument rc_eff requires input arguments catm, ra and rb'
+         !END IF
+         ! compute rc_eff :
+         !CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
+         !ENDIF
+      ENDIF
+
+    END SUBROUTINE drydepos_gas_depac
+
+
+
+!-------------------------------------------------------------------
+! rc_special: compute total canopy resistance in special CASEs
+!-------------------------------------------------------------------
+    SUBROUTINE rc_special(icmp,compnam,lu,t,nwet,rc_tot,ready,ccomp_tot)
+
+      INTEGER(iwp)         , INTENT(in)  :: icmp            ! component index
+      CHARACTER(len=*)     , INTENT(in)  :: compnam         ! component name
+      INTEGER(iwp)         , INTENT(in)  :: lu              ! land use type, lu = 1,...,nlu
+      REAL(kind=wp)        , INTENT(in)  :: t               ! temperature (C)
+      ! = 1 -> constant Rc
+      ! = 2 -> temperature dependent Rc
+      INTEGER(iwp)         , INTENT(in)  :: nwet            ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+      REAL(kind=wp)        , INTENT(out) :: rc_tot          ! total canopy resistance Rc (s/m)
+      LOGICAL              , INTENT(out) :: ready           ! Rc has been set
+      REAL(kind=wp)        , INTENT(out) :: ccomp_tot       ! total compensation point (ug/m3)
+
+      ! rc_tot is not yet set:
+      ready = .false.
+
+      ! Default compensation point in special CASEs = 0:
+      ccomp_tot = 0.0
+
+      SELECT CASE(trim(compnam))
+      CASE('HNO3','N2O5','NO3','H2O2')
+         ! No separate resistances for HNO3; just one total canopy resistance:
+         IF (t .lt. -5.0 .and. nwet .eq. 9) then
+            ! T < 5 C and snow:
+            rc_tot = 50.
+         ELSE
+            ! all other circumstances:
+            rc_tot = 10.0
+         ENDIF
+         ready = .true.
+
+      CASE('NO','CO')
+         IF (lu .eq. ilu_water_sea .or. lu .eq. ilu_water_inland) then       ! water
+            rc_tot = 2000.
+            ready = .true.
+         ELSEIF (nwet .eq. 1) then ! wet
+            rc_tot = 2000.
+            ready = .true.
+         ENDIF
+      CASE('NO2','O3','SO2','NH3')
+         ! snow surface:
+         IF (nwet.eq.9) then
+            !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
+            ready = .true.
+         ENDIF
+      CASE default
+         message_string = 'Component '// trim(compnam) // ' not supported'
+         CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 )
+      END SELECT
+
+    END SUBROUTINE rc_special
+
+
+
+!-------------------------------------------------------------------
+! rc_gw: compute external conductance
+!-------------------------------------------------------------------
+    SUBROUTINE rc_gw( compnam, iratns,t,rh,nwet, SAI_present, sai, gw )
+
+      ! Input/output variables:
+      CHARACTER(len=*)     , INTENT(in)  :: compnam ! component name
+      INTEGER(iwp)         , INTENT(in)  :: iratns  ! index for NH3/SO2 ratio;
+      ! iratns = 1: low NH3/SO2
+      ! iratns = 2: high NH3/SO2
+      ! iratns = 3: very low NH3/SO2
+      REAL(kind=wp)         , INTENT(in)  :: t       ! temperature (C)
+      REAL(kind=wp)         , INTENT(in)  :: rh      ! relative humidity (%)
+      INTEGER(iwp)          , INTENT(in)  :: nwet    ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+      LOGICAL               , INTENT(in)  :: SAI_present
+      REAL(kind=wp)         , INTENT(in)  :: sai     ! one-sided leaf area index (-)
+      REAL(kind=wp)         , INTENT(out) :: gw      ! external leaf conductance (m/s)
+
+      SELECT CASE(trim(compnam))
+
+         !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3
+
+      CASE('NO2')
+         CALL rw_constant(2000.0_wp,SAI_present,gw)
+
+      CASE('NO','CO')
+         CALL rw_constant(-9999.0_wp,SAI_present,gw)  ! see Erisman et al, 1994 section 3.2.3
+
+      CASE('O3')
+         CALL rw_constant(2500.0_wp,SAI_present,gw)
+
+      CASE('SO2')
+         CALL rw_so2( t, nwet, rh, iratns, SAI_present, gw )
+
+      CASE('NH3')
+         CALL rw_nh3_sutton(t,rh,SAI_present,gw)
+
+         ! conversion from leaf resistance to canopy resistance by multiplying with SAI:
+         Gw = sai*gw
+
+      CASE default
+         message_string = 'Component '// trim(compnam) // ' not supported'
+         CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 )
+      END SELECT
+
+    END SUBROUTINE rc_gw
+
+
+
+!-------------------------------------------------------------------
+! rw_so2: compute external leaf conductance for SO2
+!-------------------------------------------------------------------
+    SUBROUTINE rw_so2( t, nwet, rh, iratns, SAI_present, gw )
+
+      ! Input/output variables:
+      REAL(kind=wp)   , INTENT(in)  :: t      ! temperature (C)
+      INTEGER(iwp)    , INTENT(in)  :: nwet   ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+      REAL(kind=wp)   , INTENT(in)  :: rh     ! relative humidity (%)
+      INTEGER(iwp)    , INTENT(in)  :: iratns ! index for NH3/SO2 ratio;
+      ! iratns = 1: low NH3/SO2
+      ! iratns = 2: high NH3/SO2
+      ! iratns = 3: very low NH3/SO2
+      LOGICAL, INTENT(in)  :: SAI_present
+      REAL(kind=wp)   , INTENT(out) :: gw     ! external leaf conductance (m/s)
+
+      ! Variables from module:
+      ! SAI_present: vegetation is present
+
+      ! Local variables:
+      REAL(kind=wp)                 :: rw     ! external leaf resistance (s/m)
+
+      ! Check IF vegetation present:
+      IF (SAI_present) then
+
+         IF (nwet .eq. 0) then
+            !--------------------------
+            ! dry surface
+            !--------------------------
+            ! T > -1 C
+            IF (t .gt. -1.0) then
+               IF (rh .lt. 81.3) then
+                  rw = 25000*exp(-0.0693*rh)
+               ELSE
+                  rw = 0.58e12*exp(-0.278*rh) + 10.
+               ENDIF
+            ELSE
+               ! -5 C < T <= -1 C
+               IF (t .gt. -5.0) then
+                  rw=200
+               ELSE
+                  ! T <= -5 C
+                  rw=500
+               ENDIF
+            ENDIF
+         ELSE
+            !--------------------------
+            ! wet surface
+            !--------------------------
+            rw = 10. !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10
+         ENDIF
+
+         ! very low NH3/SO2 ratio:
+         IF (iratns == iratns_very_low ) rw = rw+50.
+
+         ! Conductance:
+         gw = 1./rw
+      ELSE
+         ! no vegetation:
+         gw = 0.0
+      ENDIF
+
+    END SUBROUTINE rw_so2
+
+
+
+!-------------------------------------------------------------------
+! rw_nh3_sutton: compute external leaf conductance for NH3,
+!                  following Sutton & Fowler, 1993
+!-------------------------------------------------------------------
+    SUBROUTINE rw_nh3_sutton(tsurf,rh,SAI_present,gw)
+
+      ! Input/output variables:
+      REAL(kind=wp)   , INTENT(in)  :: tsurf     ! surface temperature (C)
+      REAL(kind=wp)   , INTENT(in)  :: rh        ! relative humidity (%)
+      LOGICAL, INTENT(in)  :: SAI_present
+      REAL(kind=wp)   , INTENT(out) :: gw        ! external leaf conductance (m/s)
+
+      ! Variables from module:
+      ! SAI_present: vegetation is present
+
+      ! Local variables:
+      REAL(kind=wp)                 :: rw                ! external leaf resistance (s/m)
+      REAL(kind=wp)                 :: sai_grass_haarweg ! surface area index at experimental site Haarweg
+
+      ! Fix SAI_grass at value valid for Haarweg data for which gamma_w parametrization is derived
+      sai_grass_haarweg = 3.5
+
+      !                  100 - rh
+      !  rw = 2.0 * exp(----------)
+      !                    12
+
+      IF (SAI_present) then
+
+         ! External resistance according to Sutton & Fowler, 1993
+         rw = 2.0 * exp((100.0 - rh)/12.0)
+         rw = sai_grass_haarweg * rw
+
+         ! Frozen soil (from Depac v1):
+         IF (tsurf .lt. 0) rw = 200
+
+         ! Conductance:
+         gw = 1./rw
+      ELSE
+         ! no vegetation:
+         gw = 0.0
+      ENDIF
+
+    END SUBROUTINE rw_nh3_sutton
+
+
+
+!-------------------------------------------------------------------
+! rw_constant: compute constant external leaf conductance
+!-------------------------------------------------------------------
+    SUBROUTINE rw_constant(rw_val,SAI_present,gw)
+
+      ! Input/output variables:
+      REAL(kind=wp)   , INTENT(in)  :: rw_val      ! constant value of Rw
+      LOGICAL         , INTENT(in)  :: SAI_present
+      REAL(kind=wp)   , INTENT(out) :: gw          ! wernal leaf conductance (m/s)
+
+      ! Variables from module:
+      ! SAI_present: vegetation is present
+
+      ! Compute conductance:
+      IF (SAI_present .and. .not.missing(rw_val)) then
+         gw = 1./rw_val
+      ELSE
+         gw = 0.
+      ENDIF
+
+    END SUBROUTINE rw_constant
+
+
+
+!-------------------------------------------------------------------
+! rc_gstom: compute stomatal conductance
+!-------------------------------------------------------------------
+    SUBROUTINE rc_gstom( icmp, compnam, lu, LAI_present, lai, glrad, sinphi, t, rh, diffc, &
+         gstom, &
+         p)
+
+      ! input/output variables:
+      INTEGER(iwp),           INTENT(in) :: icmp           ! component index
+      CHARACTER(len=*),       INTENT(in) :: compnam        ! component name
+      INTEGER(iwp),           INTENT(in) :: lu             ! land use type , lu = 1,...,nlu
+      LOGICAL,                INTENT(in) :: LAI_present
+      REAL(kind=wp),          INTENT(in) :: lai            ! one-sided leaf area index
+      REAL(kind=wp),          INTENT(in) :: glrad          ! global radiation (W/m2)
+      REAL(kind=wp),          INTENT(in) :: sinphi         ! sin of solar elevation angle
+      REAL(kind=wp),          INTENT(in) :: t              ! temperature (C)
+      REAL(kind=wp),          INTENT(in) :: rh             ! relative humidity (%)
+      REAL(kind=wp),          INTENT(in) :: diffc              ! diffusion coefficient of the gas involved
+      REAL(kind=wp),          INTENT(out):: gstom              ! stomatal conductance (m/s)
+      REAL(kind=wp), OPTIONAL,INTENT(in) :: p              ! pressure (Pa)
+
+
+      ! Local variables
+      REAL(kind=wp)                      :: vpd            ! vapour pressure deficit (kPa)
+
+      REAL(kind=wp), PARAMETER           :: dO3 = 0.13e-4  ! diffusion coefficient of ozon (m2/s)
+
+
+      SELECT CASE(trim(compnam))
+
+         !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3
+
+      CASE('NO','CO')
+         ! for NO stomatal uptake is neglected:
+         gstom = 0.0
+
+      CASE('NO2','O3','SO2','NH3')
+
+         ! IF vegetation present:
+         IF (LAI_present) then
+
+            IF (glrad .gt. 0.0) then
+               CALL rc_get_vpd(t,rh,vpd)
+               CALL rc_gstom_emb( lu, glrad, t, vpd, LAI_present, lai, sinphi, gstom, p )
+               gstom = gstom*diffc/dO3       ! Gstom of Emberson is derived for ozone 
+            ELSE
+               gstom = 0.0
+            ENDIF
+         ELSE
+            ! no vegetation; zero conductance (infinite resistance):
+            gstom = 0.0
+         ENDIF
+
+      CASE default
+         message_string = 'Component '// trim(compnam) // ' not supported'
+         CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 )
+      END SELECT
+
+    END SUBROUTINE rc_gstom
+
+
+
+!-------------------------------------------------------------------
+! rc_gstom_emb: stomatal conductance according to Emberson
+!-------------------------------------------------------------------
+    SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, LAI_present, lai, sinp, &
+         Gsto, p )
+      ! History
+      !   Original code from Lotos-Euros, TNO, M. Schaap
+      !   2009-08, M.C. van Zanten, Rivm
+      !     Updated and extended.
+      !   2009-09, Arjo Segers, TNO
+      !     Limitted temperature influence to range to avoid
+      !     floating point exceptions.
+      !
+      ! Method
+      !
+      !   Code based on Emberson et al, 2000, Env. Poll., 403-413
+      !   Notation conform Unified EMEP Model Description Part 1, ch 8
+      !
+      !   In the calculation of F_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
+      !   parametrizations of Norman 1982 are applied
+      !
+      !   F_phen and F_SWP are set to 1
+      !
+      !   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
+      !   DEPAC                     Emberson
+      !     1 = grass                 GR = grassland
+      !     2 = arable land           TC = temperate crops ( LAI according to RC = rootcrops)
+      !     3 = permanent crops       TC = temperate crops ( LAI according to RC = rootcrops)
+      !     4 = coniferous forest     CF = tempareate/boREAL(kind=wp) coniferous forest
+      !     5 = deciduous forest      DF = temperate/boREAL(kind=wp) deciduous forest
+      !     6 = water                 W  = water
+      !     7 = urban                 U  = urban
+      !     8 = other                 GR = grassland
+      !     9 = desert                DE = desert
+      !
+
+      ! Emberson specific declarations
+
+      ! Input/output variables:
+      INTEGER(iwp),               INTENT(in) :: lu                ! land use type, lu = 1,...,nlu
+      REAL(kind=wp),              INTENT(in) :: glrad             ! global radiation (W/m2)
+      REAL(kind=wp),              INTENT(in) :: T                 ! temperature (C)
+      REAL(kind=wp),              INTENT(in) :: vpd               ! vapour pressure deficit (kPa)
+      LOGICAL,                    INTENT(in) :: LAI_present
+      REAL(kind=wp),              INTENT(in) :: lai               ! one-sided leaf area index
+      REAL(kind=wp),              INTENT(in) :: sinp              ! sin of solar elevation angle
+      REAL(kind=wp),              INTENT(out):: Gsto              ! stomatal conductance (m/s)
+      REAL(kind=wp), OPTIONAL,    INTENT(in) :: p                 ! pressure (Pa)
+
+      ! Local variables:
+      REAL(kind=wp)                          ::  F_light, F_phen, F_temp, F_vpd, F_swp, bt, F_env
+      REAL(kind=wp)                          ::  pardir, pardiff, parshade, parsun, LAIsun, LAIshade, sinphi
+      REAL(kind=wp)                          ::  pres
+      REAL(kind=wp), PARAMETER               ::  p_sealevel = 1.01325e5    ! Pa
+
+
+      ! Check whether vegetation is present:
+      IF (LAI_present) then
+
+         ! calculation of correction factors for stomatal conductance
+         IF (sinp .le. 0.0) then  
+            sinphi = 0.0001
+         ELSE
+            sinphi = sinp
+         END IF
+
+         ! ratio between actual and sea-level pressure is used
+         ! to correct for height in the computation of par ;
+         ! should not exceed sea-level pressure therefore ...
+         IF ( present(p) ) then
+            pres = min( p, p_sealevel )
+         ELSE 
+            pres = p_sealevel
+         ENDIF
+
+         ! Direct and diffuse par, Photoactive (=visible) radiation:
+         CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff )
+
+         ! par for shaded leaves (canopy averaged):
+         parshade = pardiff*exp(-0.5*LAI**0.7)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Norman, 1982
+         IF (glrad .gt. 200 .and. LAI .gt. 2.5) then
+            parshade = pardiff*exp(-0.5*LAI**0.8)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Zhang et al., 2001
+         END IF
+
+         ! par for sunlit leaves (canopy averaged):
+         ! alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
+         parsun = pardir*0.5/sinphi + parshade ! Norman, 1982
+         IF (glrad .gt. 200 .and. LAI .gt. 2.5) then
+            parsun = pardir**0.8*0.5/sinphi + parshade ! Zhang et al., 2001
+         END IF
+
+         ! leaf area index for sunlit and shaded leaves:
+         IF (sinphi .gt. 0) then
+            LAIsun = 2*sinphi*(1-exp(-0.5*LAI/sinphi ))
+            LAIshade = LAI -LAIsun
+         ELSE
+            LAIsun = 0
+            LAIshade = LAI
+         END IF
+
+         F_light = (LAIsun*(1 - exp(-1.*alpha(lu)*parsun)) + LAIshade*(1 - exp(-1.*alpha(lu)*parshade)))/LAI 
+
+         F_light = max(F_light,F_min(lu))
+
+         !  temperature influence; only non-zero within range [Tmin,Tmax]:
+         IF ( (Tmin(lu) < T) .and. (T < Tmax(lu)) ) then
+            BT = (Tmax(lu)-Topt(lu))/(Topt(lu)-Tmin(lu))
+            F_temp = ((T-Tmin(lu))/(Topt(lu)-Tmin(lu))) * ((Tmax(lu)-T)/(Tmax(lu)-Topt(lu)))**BT
+         ELSE
+            F_temp = 0.0
+         END IF
+         F_temp = max(F_temp,F_min(lu))
+
+         ! vapour pressure deficit influence
+         F_vpd = min(1.,((1.-F_min(lu))*(vpd_min(lu)-vpd)/(vpd_min(lu)-vpd_max(lu)) + F_min(lu) ))
+         F_vpd = max(F_vpd,F_min(lu))
+
+         F_swp = 1.
+
+         ! influence of phenology on stom. conductance
+         ! ignored for now in DEPAC since influence of F_phen on lu classes in use is negligible.
+         ! When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
+         F_phen = 1.
+
+         ! evaluate total stomatal conductance
+         F_env = F_temp*F_vpd*F_swp
+         F_env = max(F_env,F_min(lu))
+         gsto = G_max(lu) * F_light * F_phen * F_env
+
+         ! gstom expressed per m2 leafarea;
+         ! this is converted with LAI to m2 surface.
+         Gsto = lai*gsto    ! in m/s
+
+      ELSE
+         GSto = 0.0
+      ENDIF
+
+    END SUBROUTINE rc_gstom_emb
+
+
+
+!-------------------------------------------------------------------
+! par_dir_diff
+!-------------------------------------------------------------------
+    SUBROUTINE par_dir_diff(glrad,sinphi,pres,pres_0,par_dir,par_diff)
+      !
+      !     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and
+      !     diffuse, visible and near-infrared components. Agric. Forest Meteorol.
+      !     34, 205-213.
+      !
+      !     From a SUBROUTINE obtained from Leiming Zhang (via Willem Asman),
+      !     MeteoroLOGICAL Service of Canada
+      !     e-mail: leiming.zhang@ec.gc.ca
+      !
+      !     Leiming uses solar irradiance. This should be equal to global radiation and
+      !     Willem Asman set it to global radiation
+
+      REAL(kind=wp),    INTENT(in)  :: glrad             ! global radiation (W m-2)
+      REAL(kind=wp),    INTENT(in)  :: sinphi            ! sine of the solar elevation
+      REAL(kind=wp),    INTENT(in)  :: pres              ! actual pressure (to correct for height) (Pa)
+      REAL(kind=wp),    INTENT(in)  :: pres_0            ! pressure at sea level (Pa)
+      REAL(kind=wp),    INTENT(out) :: par_dir           ! PAR direct : visible (photoactive) direct beam radiation (W m-2)
+      REAL(kind=wp),    INTENT(out) :: par_diff          ! PAR diffuse: visible (photoactive) diffuse radiation (W m-2)
+
+      !     fn              = near-infrared direct beam fraction (DIMENSIONless)
+      !     fv              = PAR direct beam fraction (DIMENSIONless)
+      !     ratio           = ratio measured to potential solar radiation (DIMENSIONless)
+      !     rdm             = potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky
+      !     rdn             = potential diffuse near-infrared radiation (W m-2)
+      !     rdu             = visible (PAR) direct beam radiation (W m-2)
+      !     rdv             = potential visible (PAR) diffuse radiation (W m-2)
+      !     rn              = near-infrared radiation (W m-2)
+      !     rv              = visible radiation (W m-2)
+      !     ww              = water absorption in the near infrared for 10 mm of precipitable water
+
+      REAL(kind=wp) :: rdu,rdv,ww,rdm,rdn,rv,rn,ratio,sv,fv
+
+      !     Calculate visible (PAR) direct beam radiation
+      !     600 W m-2 represents average amount of PAR (400-700 nm wavelength)
+      !     at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
+      rdu=600.*exp(-0.185*(pres/pres_0)/sinphi)*sinphi
+
+      !     Calculate potential visible diffuse radiation
+      rdv=0.4*(600.- rdu)*sinphi
+
+      !     Calculate the water absorption in the-near infrared
+      ww=1320*10**( -1.195+0.4459*log10(1./sinphi)-0.0345*(log10(1./sinphi))**2 )
+
+      !     Calculate potential direct beam near-infrared radiation
+      rdm=(720.*exp(-0.06*(pres/pres_0)/sinphi)-ww)*sinphi     !720 = solar constant - 600
+
+      !     Calculate potential diffuse near-infrared radiation
+      rdn=0.6*(720-rdm-ww)*sinphi
+
+      ! Compute visible and near-infrared radiation
+      rv=max(0.1,rdu+rdv)
+      rn=max(0.01,rdm+rdn)
+
+      ! Compute ratio between input global radiation and total radiation computed here
+      ratio=min(0.9,glrad/(rv+rn))
+
+      !     Calculate total visible radiation
+      sv=ratio*rv
+
+      !     Calculate fraction of PAR in the direct beam
+      fv=min(0.99, (0.9-ratio)/0.7)            ! help variable
+      fv=max(0.01,rdu/rv*(1.0-fv**0.6667))     ! fraction of PAR in the direct beam
+
+      ! Compute direct and diffuse parts of PAR
+      par_dir=fv*sv
+      par_diff=sv-par_dir
+
+    END SUBROUTINE par_dir_diff
+
+!-------------------------------------------------------------------
+! rc_get_vpd: get vapour pressure deficit (kPa)
+!-------------------------------------------------------------------
+    SUBROUTINE rc_get_vpd(temp, relh,vpd)
+
+      IMPLICIT NONE
+
+      ! Input/output variables:
+      REAL(kind=wp)    , INTENT(in)  :: temp   ! temperature (C)
+      REAL(kind=wp)    , INTENT(in)  :: relh   ! relative humidity (%)
+      REAL(kind=wp)    , INTENT(out) :: vpd    ! vapour pressure deficit (kPa)
+
+      ! Local variables:
+      REAL(kind=wp) :: esat
+
+      ! fit PARAMETERs:
+      REAL(kind=wp), PARAMETER :: a1 = 6.113718e-1
+      REAL(kind=wp), PARAMETER :: a2 = 4.43839e-2
+      REAL(kind=wp), PARAMETER :: a3 = 1.39817e-3
+      REAL(kind=wp), PARAMETER :: a4 = 2.9295e-5
+      REAL(kind=wp), PARAMETER :: a5 = 2.16e-7
+      REAL(kind=wp), PARAMETER :: a6 = 3.0e-9
+
+      ! esat is saturation vapour pressure (kPa) at temp(C) following monteith(1973)
+      esat = a1 + a2*temp + a3*temp**2 + a4*temp**3 + a5*temp**4 + a6*temp**5
+      vpd  = esat * (1-relh/100)
+
+    END SUBROUTINE rc_get_vpd
+
+
+
+!-------------------------------------------------------------------
+! rc_gsoil_eff: compute effective soil conductance
+!-------------------------------------------------------------------
+    SUBROUTINE rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff)
+
+      ! Input/output variables:
+      INTEGER(iwp)    , INTENT(in)  :: icmp           ! component index
+      INTEGER(iwp)    , INTENT(in)  :: lu             ! land use type, lu = 1,...,nlu
+      REAL(kind=wp)   , INTENT(in)  :: sai            ! surface area index
+      REAL(kind=wp)   , INTENT(in)  :: ust            ! friction velocity (m/s)
+      INTEGER(iwp)    , INTENT(in)  :: nwet           ! index for wetness
+      ! nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+      ! N.B. this routine cannot be CALLed with nwet = 9,
+      ! nwet = 9 should be handled outside this routine.
+      REAL(kind=wp)   , INTENT(in)  :: t              ! temperature (C)
+      REAL(kind=wp)   , INTENT(out) :: gsoil_eff      ! effective soil conductance (m/s)
+
+      ! local variables:
+      REAL(kind=wp)                 :: rinc           ! in canopy resistance  (s/m)
+      REAL(kind=wp)                 :: rsoil_eff      ! effective soil resistance (s/m)
+
+
+
+      ! Soil resistance (numbers matched with lu_classes and component numbers)
+      REAL(kind=wp), PARAMETER     :: rsoil(nlu_dep,ncmp) = reshape( (/ &
+                                ! grs    ara    crp    cnf    dec    wat    urb    oth    des    ice    sav    trf    wai    med    sem
+           1000.,  200.,  200.,  200.,  200., 2000.,  400., 1000., 2000., 2000., 1000.,  200., 2000.,  200.,  400., & ! O3
+           1000., 1000., 1000., 1000., 1000.,   10., 1000., 1000., 1000.,  500., 1000., 1000.,   10., 1000., 1000., & ! SO2
+           1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., & ! NO2
+           -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! NO
+           100.,  100.,  100.,  100.,  100.,   10.,  100.,  100.,  100., 1000.,  100.,  100.,   10.,  100.,  100.,  & ! NH3
+           -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! CO 
+           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! NO3 
+           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! HNO3
+           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! N2O5
+           -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),& ! H2O2   
+           (/nlu_dep,ncmp/) )
+
+      !  o3       so2       no2        no       nh3        co    no3   hno3   n2o5   h2o2
+
+
+      REAL(kind=wp), PARAMETER     :: rsoil_wet(ncmp)    = (/2000.,      10.,    2000.,    -999.,      10.,    -999., -999., -999., -999., -999./)
+      REAL(kind=wp), PARAMETER     :: rsoil_frozen(ncmp) = (/2000.,     500.,    2000.,    -999.,    1000.,    -999., -999., -999., -999., -999./)
+
+
+
+      ! Compute in canopy (in crop) resistance:
+      CALL rc_rinc(lu,sai,ust,rinc)
+
+      ! Check for missing deposition path:
+      IF (missing(rinc)) then
+         rsoil_eff = -9999.
+      ELSE
+         ! Frozen soil (temperature below 0):
+         IF (t .lt. 0.0) then
+            IF (missing(rsoil_frozen(icmp))) then
+               rsoil_eff = -9999.
+            ELSE
+               rsoil_eff = rsoil_frozen(icmp) + rinc
+            ENDIF
+         ELSE
+            ! Non-frozen soil; dry:
+            IF (nwet .eq. 0) then
+               IF (missing(rsoil(lu,icmp))) then
+                  rsoil_eff = -9999.
+               ELSE
+                  rsoil_eff = rsoil(lu,icmp) + rinc
+               ENDIF
+
+               ! Non-frozen soil; wet:
+            ELSEIF (nwet .eq. 1) then
+               IF (missing(rsoil_wet(icmp))) then
+                  rsoil_eff = -9999.
+               ELSE
+                  rsoil_eff = rsoil_wet(icmp) + rinc
+               ENDIF
+            ELSE
+               message_string = 'nwet can only be 0 or 1'
+               CALL message( 'rc_gsoil_eff', 'CM0460', 1, 2, 0, 6, 0 )
+            ENDIF
+         ENDIF
+      ENDIF
+
+      ! Compute conductance:
+      IF (rsoil_eff .gt. 0.0) then
+         gsoil_eff = 1./rsoil_eff
+      ELSE
+         gsoil_eff = 0.0
+      ENDIF
+
+    END SUBROUTINE rc_gsoil_eff
+
+
+
+!-------------------------------------------------------------------
+! rc_rinc: compute in canopy (or in crop) resistance
+! van Pul and Jacobs, 1993, BLM
+!-------------------------------------------------------------------
+    SUBROUTINE rc_rinc(lu,sai,ust,rinc)
+
+
+      ! Input/output variables:
+      INTEGER(iwp)    , INTENT(in)  :: lu          ! land use class, lu = 1, ..., nlu
+      REAL(kind=wp)   , INTENT(in)  :: sai         ! surface area index
+      REAL(kind=wp)   , INTENT(in)  :: ust         ! friction velocity (m/s)
+      REAL(kind=wp)   , INTENT(out) :: rinc        ! in canopy resistance (s/m)
+
+
+      ! b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 IF not applicable)
+      ! h = vegetation height (m)                             grass arabl  crops conIF decid   water  urba   othr  desr   ice  sav  trf   wai   med semi
+      REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: b = (/ -999,  14,    14,   14,   14,    -999,  -999, -999, -999, -999, -999, 14, -999, 14, 14 /)
+      REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: h = (/ -999,  1,     1,    20,   20,    -999,  -999, -999, -999, -999, -999, 20,  -999, 1 ,  1 /)
+
+      ! Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
+      IF (b(lu) .gt. 0.0) then
+
+         ! Check for u* > 0 (otherwise denominator = 0):
+         IF (ust .gt. 0.0) then
+            rinc = b(lu)*h(lu)*sai/ust
+         ELSE
+            rinc = 1000.0
+         ENDIF
+      ELSE
+         IF (lu .eq. ilu_grass .or. lu .eq. ilu_other ) then
+            rinc = -999. ! no deposition path for grass, other, and semi-natural
+         ELSE
+            rinc = 0.0   ! no in-canopy resistance
+         ENDIF
+      ENDIF
+
+    END SUBROUTINE rc_rinc
+
+
+
+!-------------------------------------------------------------------
+! rc_rctot: compute total canopy (or surface) resistance Rc
+!-------------------------------------------------------------------
+    SUBROUTINE rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot)
+
+      ! Input/output variables:
+      REAL(kind=wp), INTENT(in)  :: gstom        ! stomatal conductance (s/m)
+      REAL(kind=wp), INTENT(in)  :: gsoil_eff    ! effective soil conductance (s/m)
+      REAL(kind=wp), INTENT(in)  :: gw           ! external leaf conductance (s/m)
+      REAL(kind=wp), INTENT(out) :: gc_tot       ! total canopy conductance (m/s)
+      REAL(kind=wp), INTENT(out) :: rc_tot       ! total canopy resistance Rc (s/m)
+
+      ! Total conductance:
+      gc_tot = gstom + gsoil_eff + gw
+
+      ! Total resistance (note: gw can be negative, but no total emission allowed here):
+      IF (gc_tot .le. 0.0 .or. gw .lt. 0.0) then
+         rc_tot = -9999.
+      ELSE
+         rc_tot = 1./gc_tot
+      ENDIF
+
+    END SUBROUTINE rc_rctot
+
+
+
+    !-------------------------------------------------------------------
+    ! rc_comp_point_rc_eff: calculate the effective resistance Rc
+    ! based on one or more compensation points
+    !-------------------------------------------------------------------
+    !
+    ! NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173.
+    ! Sutton 1998 AE 473-480)
+    !
+    ! Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine.
+    ! FS 2009-01-29: variable names made consistent with DEPAC 
+    ! FS 2009-03-04: use total compensation point
+    !
+    ! C: with total compensation point   ! D: approximation of C
+    !                                    !    with classical approach
+    !  zr --------- Catm                 !  zr --------- Catm    
+    !         |                          !         |       
+    !         Ra                         !         Ra      
+    !         |                          !         |       
+    !         Rb                         !         Rb      
+    !         |                          !         |       
+    !  z0 --------- Cc                   !  z0 --------- Cc
+    !         |                          !         |              
+    !        Rc                          !        Rc_eff          
+    !         |                          !         |              
+    !     --------- Ccomp_tot            !     --------- C=0
+    !
+    !
+    ! The effective Rc is defined such that instead of using
+    !
+    !   F = -vd*[Catm - Ccomp_tot]                                    (1)
+    !
+    ! we can use the 'normal' flux formula
+    !
+    !   F = -vd'*Catm,                                                (2)
+    !
+    ! with vd' = 1/(Ra + Rb + Rc')                                    (3)
+    !
+    ! and Rc' the effective Rc (rc_eff).
+    !                                                (Catm - Ccomp_tot)
+    ! vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------
+    !                                                      Catm
+    !
+    !                                        (Catm - Ccomp_tot)
+    ! 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------
+    !                                              Catm
+    !
+    !                                          Catm
+    ! (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------
+    !                                     (Catm - Ccomp_tot)
+    !
+    !                              Catm
+    ! Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb
+    !                        (Catm - Ccomp_tot)
+    !
+    !                        Catm                           Catm
+    ! Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------
+    !                  (Catm - Ccomp_tot)              (Catm - Ccomp_tot)
+    !
+    ! Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot)
+    !
+    ! This is not the most efficient way to DO this;
+    ! in the current LE version, (1) is used directly in the CALLing routine
+    ! and this routine is not used anymore.
+    !
+    ! -------------------------------------------------------------------------------------------
+    ! In fact it is the question IF this correct; note the difference in differential equations
+    !
+    !   dCatm                              !   dCatm
+    ! H ----- = F = -vd'*Catm,             ! H ----- = F = -vd*[Catm - Ccomp_tot],
+    !    dt                                !    dt
+    !
+    ! where in the left colum vd' is a function of Catm and not a constant anymore.
+    !
+    ! Another problem is that IF Catm -> 0, we end up with an infinite value of the exchange
+    ! velocity vd.
+    ! -------------------------------------------------------------------------------------------
+
+    SUBROUTINE rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
+
+      IMPLICIT NONE
+
+      ! Input/output variables:
+      REAL(kind=wp), INTENT(in)  :: ccomp_tot  ! total compensation point (weighed average of separate compensation points) (ug/m3)
+      REAL(kind=wp), INTENT(in)  :: catm       ! atmospheric concentration (ug/m3)
+      REAL(kind=wp), INTENT(in)  :: ra         ! aerodynamic resistance (s/m)
+      REAL(kind=wp), INTENT(in)  :: rb         ! boundary layer resistance (s/m)
+      REAL(kind=wp), INTENT(in)  :: rc_tot     ! total canopy resistance (s/m)
+      REAL(kind=wp), INTENT(out) :: rc_eff     ! effective total canopy resistance (s/m)
+
+      ! Compute effective resistance:
+      IF ( ccomp_tot == 0.0 ) then
+         ! trace with no compensiation point ( or compensation point equal to zero)
+         rc_eff = rc_tot
+
+      ELSE IF ( ccomp_tot > 0 .and. ( abs( catm - ccomp_tot ) < 1.e-8 ) ) then
+         ! surface concentration (almoast) equal to atmospheric concentration
+         ! no exchange between surface and atmosphere, infinite RC --> vd=0
+         rc_eff = 9999999999.
+
+      ELSE IF ( ccomp_tot > 0 ) then
+         ! compensation point available, calculate effective restisance
+         rc_eff = ((ra + rb)*ccomp_tot + rc_tot*catm)/(catm-ccomp_tot)
+
+      ELSE
+         rc_eff = -999.
+         message_string = 'This should not be possible, check ccomp_tot'
+         CALL message( 'rc_comp_point_rc_eff', 'CM0461', 1, 2, 0, 6, 0 )
+      ENDIF
+
+      return
+    END SUBROUTINE rc_comp_point_rc_eff
+
+
+
+!-------------------------------------------------------------------
+! missing: check for data that correspond with a missing deposition path
+!          this data is represented by -999
+!-------------------------------------------------------------------
+
+    LOGICAL function missing(x)
+
+      REAL(kind=wp), INTENT(in) :: x
+
+      ! bandwidth for checking (in)equalities of floats
+      REAL(kind=wp), PARAMETER :: EPS = 1.0e-5
+
+      missing = (abs(x + 999.) .le. EPS)
+
+    END function missing
+
+
+
+    ELEMENTAL FUNCTION sedimentation_velocity(rhopart, partsize, slipcor, visc) RESULT (vs)
+
+
+      IMPLICIT NONE
+
+      ! --- in/out ---------------------------------
+
+      REAL(kind=wp), INTENT(in)        ::  rhopart    ! kg/m3
+      REAL(kind=wp), INTENT(in)        ::  partsize   ! m
+      REAL(kind=wp), INTENT(in)        ::  slipcor    ! m
+      REAL(kind=wp), INTENT(in)        ::  visc       ! viscosity
+      REAL(kind=wp)                    ::  vs
+
+      ! acceleration of gravity:
+      REAL(kind=wp), PARAMETER         ::  grav = 9.80665    ! m/s2
+
+      ! --- begin ----------------------------------
+
+      !viscosity & sedimentation velocity
+      vs = rhopart * (partsize**2) * grav * slipcor / (18*visc)
+
+    END FUNCTION sedimentation_velocity
+
+!------------------------------------------------------------------------
+! Boundary-layer deposition resistance following Zhang (2001)
+!------------------------------------------------------------------------
+    SUBROUTINE drydepo_aero_zhang_vd( vd, Rs, vs1, partsize, slipcor, &
+         nwet, tsurf, dens1, viscos1, &
+         luc, ftop_lu, ustar_lu)
+
+
+      IMPLICIT NONE  
+
+      ! --- in/out ---------------------------------
+
+      REAL(kind=wp), INTENT(out)        ::  vd          ! deposition velocity (m/s)
+      REAL(kind=wp), INTENT(out)        ::  Rs          ! Sedimentaion resistance (s/m)
+      REAL(kind=wp), INTENT(in)         ::  vs1         ! sedimentation velocity in lowest layer
+      REAL(kind=wp), INTENT(in)         ::  partsize    ! particle diameter (m)
+      REAL(kind=wp), INTENT(in)         ::  slipcor     ! slip correction factor
+      REAL(kind=wp), INTENT(in)         ::  tsurf       ! surface temperature (K)
+      REAL(kind=wp), INTENT(in)         ::  dens1       ! air density (kg/m3) in lowest layer
+      REAL(kind=wp), INTENT(in)         ::  viscos1     ! air viscosity in lowest layer
+      REAL(kind=wp), INTENT(in)         ::  ftop_lu     ! atmospheric resistnace Ra
+      REAL(kind=wp), INTENT(in)         ::  ustar_lu    ! friction velocity u*    
+
+      INTEGER(iwp),  INTENT(in)         ::  nwet        ! 1=rain, 9=snowcover
+      INTEGER(iwp),  INTENT(in)         ::  luc         ! DEPAC LU
+
+
+      ! --- const ----------------------------------
+
+
+
+      ! acceleration of gravity:
+      REAL(kind=wp), PARAMETER         ::  grav     = 9.80665    ! m/s2
+
+      REAL(kind=wp), PARAMETER         ::  beta     = 2.0
+      REAL(kind=wp), PARAMETER         ::  epsilon0 = 3.0
+      REAL(kind=wp), PARAMETER         ::  kb       = 1.38066e-23
+      REAL(kind=wp), PARAMETER         ::  pi = 3.141592654_wp  !< PI
+
+      REAL, PARAMETER :: alfa_lu(nlu_dep) = & 
+           (/1.2,  1.2,   1.2,  1.0,  1.0,   100.0, 1.5,  1.2, 50.0, 100.0, 1.2, 1.0, 100.0, 1.2, 50.0/)   
+      !           grass arabl  crops conif decid  water  urba  othr  desr  ice   sav   trf  wai   med  sem 
+      REAL, PARAMETER :: gamma_lu(nlu_dep) = &
+           (/0.54, 0.54,  0.54, 0.56, 0.56,  0.50,  0.56, 0.54, 0.58, 0.50, 0.54, 0.56, 0.50, 0.54, 0.54/) 
+      !           grass arabl  crops conif decid  water  urba  othr  desr  ice   sav   trf   wai   med   sem   
+      REAL, PARAMETER ::A_lu(nlu_dep) = &   
+           (/3.0,  3.0,   2.0,  2.0,  7.0,  -99.,   10.0, 3.0, -99., -99.,  3.0, 7.0, -99., 2.0, -99./)
+      !           grass arabl  crops conif decid  water  urba  othr  desr  ice   sav  trf   wai  med   sem     
+
+
+      ! --- local ----------------------------------
+
+      REAL(kind=wp)              ::  kinvisc, diff_part
+      REAL(kind=wp)              ::  schmidt,stokes, Ebrown, Eimpac, Einterc, Reffic
+
+      ! --- begin ----------------------------------
+
+      ! kinetic viscosity & diffusivity
+      kinvisc = viscos1 / dens1 ! only needed at surface
+
+      diff_part = kb * tsurf * slipcor / (3*pi*viscos1*partsize)
+
+      ! Schmidt number
+      schmidt = kinvisc / diff_part
+
+      ! calculate collection efficiencie E
+      Ebrown = Schmidt**(-gamma_lu(luc)) ! Brownian diffusion
+      ! determine Stokes number, interception efficiency 
+      ! and sticking efficiency R (1 = no rebound)
+      IF ( luc == ilu_ice .or. nwet.eq.9 .or. luc == ilu_water_sea .or. luc == ilu_water_inland ) THEN
+         stokes=vs1*ustar_lu**2/(grav*kinvisc)
+         Einterc=0.0
+         Reffic=1.0
+      ELSE IF ( luc == ilu_other .or. luc == ilu_desert ) THEN !tundra of desert
+         stokes=vs1*ustar_lu**2/(grav*kinvisc)
+         Einterc=0.0
+         Reffic=exp(-Stokes**0.5)
+      ELSE
+         stokes=vs1*ustar_lu/(grav*A_lu(luc)*1.e-3)
+         Einterc=0.5*(partsize/(A_lu(luc)*1e-3))**2
+         Reffic=exp(-Stokes**0.5)
+      END IF
+      ! when surface is wet all particles DO not rebound:
+      IF(nwet.eq.1) Reffic = 1.0
+      ! determine impaction efficiency:
+      Eimpac = ( stokes / (alfa_lu(luc)+stokes) )**beta
+
+      ! sedimentation resistance:
+      Rs = 1.0 / ( epsilon0 * ustar_lu * (Ebrown+Eimpac+Einterc) * Reffic )
+
+      ! deposition velocity according to Seinfeld and Pandis (= SP, 2006; eq 19.7):
+      !
+      !            1
+      !    vd = ------------------ + vs
+      !         Ra + Rs + Ra*Rs*vs
+      !
+      !    where: Rs = Rb (in SP, 2006)
+      !
+      !
+      vd = 1.0 / ( ftop_lu + Rs + ftop_lu*Rs*vs1) + vs1
+
+
+
+    END SUBROUTINE drydepo_aero_zhang_vd
+
+
+
+    SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb )    
+
+      ! Compute quasi-laminar boundary layer resistance as a function of landuse and tracer.
+      ! 
+      ! Original EMEP formulation by (Simpson et al, 2003) is used. 
+      !
+
+
+      IMPLICIT NONE
+
+      ! --- in/out ---------------------------------
+
+      LOGICAL      , INTENT(in)      ::  is_water
+      REAL(kind=wp), INTENT(in)      ::  z0h           ! roughness length for heat
+      REAL(kind=wp), INTENT(in)      ::  ustar         ! friction velocity
+      REAL(kind=wp), INTENT(in)      ::  diffc         ! coefficient of diffusivity
+      REAL(kind=wp), INTENT(out)     ::  Rb            ! boundary layer resistance
+
+      ! --- const ----------------------------------
+
+      ! thermal diffusivity of dry air 20 C
+      REAL(kind=wp), PARAMETER       ::  thk = 0.19e-4         ! http://en.wikipedia.org/wiki/Thermal_diffusivity  (@ T=300K)
+      REAL(kind=wp), PARAMETER       ::  kappa_stab = 0.35     ! von Karman constant
+
+      !  --- begin ---------------------------------
+
+
+      ! Use Simpson et al. (2003) 
+      !IF ( is_water ) THEN
+      !Rb = 1.0_wp / (kappa_stab*max(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.01_wp,ustar))
+      !TODO: Check Rb over water calculation!!!
+      !   Rb = 1.0_wp / (kappa_stab*max(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.1_wp,ustar))
+      !ELSE
+      Rb = 5.0_wp / max(0.01_wp,ustar) * (thk/diffc)**0.67_wp
+      !END IF
+
+    END SUBROUTINE get_rb_cell
+
+
+!-----------------------------------------------------------------
+!  Compute water vapor partial pressure (e_w)
+!  given specific humidity Q [(kg water)/(kg air)].
+!
+!  Use that gas law for volume V with temperature T
+!  holds for the total mixture as well as the water part:
+!
+!    R T / V = p_air / n_air = p_water / n_water
+!
+!  thus:
+!
+!    p_water = p_air n_water / n_air
+!
+!  Use:
+!    n_air =   m_air   /        xm_air
+!            [kg air]  /  [(kg air)/(mole air)]
+!  and:
+!    n_water =  m_air * Q  /     xm_water
+!              [kg water]  /  [(kg water)/(mole water)]
+!  thus:
+!    p_water = p_air Q / (xm_water/xm_air)
+!
+
+    ELEMENTAL FUNCTION watervaporpartialpressure( Q, p ) result ( p_w )
+
+
+      ! --- in/out ---------------------------------
+
+      REAL(kind=wp), INTENT(in)      ::  Q   ! specific humidity [(kg water)/(kg air)]
+      REAL(kind=wp), INTENT(in)      ::  p   ! air pressure [Pa]
+      REAL(kind=wp)                  ::  p_w  ! water vapor partial pressure [Pa]
+
+      ! --- const ----------------------------------
+
+      ! mole mass ratio:
+      REAL(kind=wp), PARAMETER  ::  eps = xm_H2O / xm_air  !  ~ 0.622
+
+      ! --- begin ----------------------------------
+
+      ! partial pressure of water vapor:
+      p_w = p * Q / eps
+
+    END function watervaporpartialpressure
+
+
+
+!------------------------------------------------------------------    
+!  Saturation vapor pressure.
+!  From (Stull 1988, eq. 7.5.2d):
+!
+!      e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) )     [Pa]
+!
+!  where:
+!      p0 = 611.2 [Pa]   : reference pressure
+!
+!  Arguments:
+!      T  [K]  : air temperature
+!  Result:
+!      e_sat_w  [Pa]  : saturation vapor pressure
+!
+!  References:
+!      Roland B. Stull, 1988
+!      An introduction to boundary layer meteorology.
+!
+
+    ELEMENTAL FUNCTION saturationvaporpressure( T ) result( e_sat_w )
+
+
+      ! --- in/out ---------------------------------
+
+      REAL(kind=wp), INTENT(in)      ::  T        ! temperature [K]
+      REAL(kind=wp)                  ::  e_sat_w  ! saturation vapor pressure  [Pa]
+
+      ! --- const ----------------------------------
+
+      ! base pressure:
+      REAL(kind=wp), PARAMETER  ::  p0 = 611.2   ! [Pa]
+
+      ! --- begin ----------------------------------
+
+      ! saturation vapor pressure:
+      e_sat_w = p0 * exp( 17.67 * (T-273.16) / (T-29.66) )   ! [Pa]
+
+    END FUNCTION saturationvaporpressure
+
+
+
+!------------------------------------------------------------------------
+!  Relative humidity RH [%] is by definition:
+!
+!           e_w            # water vapor partial pressure
+!    Rh = -------- * 100
+!         e_sat_w          # saturation vapor pressure
+!
+!  Compute here from:
+!    Q specific humidity [(kg water)/(kg air)]
+!    T temperature [K]
+!    P pressure [Pa]
+!
+    
+    ELEMENTAL FUNCTION relativehumidity_from_specifichumidity( Q, T, p ) result( Rh )
+
+
+      ! --- in/out ---------------------------------
+
+      REAL(kind=wp), INTENT(in)      ::  Q   ! specific humidity [(kg water)/(kg air)]
+      REAL(kind=wp), INTENT(in)      ::  T   ! temperature [K]
+      REAL(kind=wp), INTENT(in)      ::  p   ! air pressure [Pa]
+      REAL(kind=wp)                  ::  Rh  ! relative humidity [%]
+
+      ! --- begin ----------------------------------
+
+      ! relative humidity:
+      Rh = watervaporpartialpressure( Q, p ) / saturationvaporpressure( T ) * 100.0
+
+    END FUNCTION relativehumidity_from_specifichumidity
+
+
+     
  END MODULE chemistry_model_mod
 

palm/trunk/SOURCE/date_and_time_mod.f90

@@ -19 +19 @@
 !
 ! Current revisions:
-! -----------------
+! ------------------
 ! 
 ! 
@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443, banzhafs:
+! Added initial hour_of_day, hour_of_year, day_of_year and month_of_year to
+! init_date_and_time
+!
+! 3298 2018-10-02 12:21:11Z kanani
 ! - Minor formatting (kanani)
 ! - Added Routines for DEFAULT mode of chemistry emissions (Russo)
@@ -61 +66 @@
 !>       already implemented changes for calculating it from date_init in 
 !>       calc_date_and_time
+!> @todo time_utc during spin-up  
 !------------------------------------------------------------------------------!
  MODULE date_and_time_mod
@@ -91 +97 @@
     INTEGER(iwp)        ::  index_mm                                !< index months of the default emission mode
     INTEGER(iwp)        ::  index_dd                                !< index days of the default emission mode
-    INTEGER(iwp)        ::  index_hh                                !< index hours of the default emission mode
+    INTEGER(iwp)        ::  index_hh                                !< index hours of the emission mode
+    INTEGER(iwp)        ::  hours_since_reference_point             !< hours of current simulation
 
     REAL(wp)            ::  time_utc                     !< current model time in UTC
@@ -102 +109 @@
     REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
     
-    CHARACTER(len=8)    ::  date_init = "31122018"                  !< Starting date of simulation: We selected this because it was a monday
+    CHARACTER(len=8)    ::  date_init = "21062017"                  !< Starting date of simulation: We selected this because it was a monday
  
     !> --- Parameters
@@ -121 +128 @@
     END INTERFACE time_default_indices
 
+    !-- Get hour index in the PRE-PROCESSED case of chemistry emissions :
+    INTERFACE time_preprocessed_indices
+       MODULE PROCEDURE time_preprocessed_indices
+    END INTERFACE time_preprocessed_indices
+
+
     !-- Calculate current date and time
     INTERFACE calc_date_and_time
@@ -128 +141 @@
 
     !-- Public Interfaces
-    PUBLIC calc_date_and_time, time_default_indices, init_date_and_time
+    PUBLIC calc_date_and_time, time_default_indices, init_date_and_time, time_preprocessed_indices
 
     !-- Public Variables
@@ -147 +160 @@
        IMPLICIT NONE
 
+       !--    Variables Definition
+       INTEGER ::  i_mon       !< Index for going through the different months
+
        IF  (day_of_year_init == 0) THEN
           ! Day of the month at starting time
@@ -159 +175 @@
        ENDIF
 
-    END SUBROUTINE init_date_and_time
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate current date and time of the simulation 
-!------------------------------------------------------------------------------!
- 
-    SUBROUTINE calc_date_and_time
-
-       IMPLICIT NONE
-
-!--    Variables Definition
-       INTEGER                          :: i_mon       !< Index for going through the different months
-
-       !> Update simulation time in seconds
-       time_update = simulated_time-coupling_start_time
-
-!--    Calculate current day of the simulated time
-       days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
-                               / 86400.0_wp ) )
-
-!--    Calculate actual UTC time                        
-       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
-       
-!sB    PRILIMINARY workaround for time_utc bug concerning time_since_reference_point:
-       time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
-
-!--    Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist. 
-
+
+       !-- Calculate initial hour of the day: the first hour of the day is from 00:00:00 to 00:59:59. 
+
+       hour_of_day = INT( FLOOR( time_utc_init/3600.0_wp ) ) + 1
+
+       !-- Calculate initial day day_of_year_init in case date_init is given or day_of_year_init is given
        IF ( day_of_year_init == 0 ) THEN
 
@@ -214 +206 @@
        ENDIF
 
+       
+       !-- Initial day of the year
+       day_of_year = day_of_year_init
+
+       !-- Initial hour of the year
+       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day
+
+       !--Initial day of the month and month of the year
+       !> --------------------------------------------------------------------------------
+       !> The first case is when date_init is not provided: we only know day_of_year_init      
+       IF ( month_of_year == 0 .AND. day_of_month == 0) THEN
+
+         
+          IF ( day_of_year .LE. 31 ) THEN
+
+             month_of_year=1
+             day_of_month=day_of_year
+
+          ELSE
+
+             DO i_mon=2,12   !january is considered in the first case
+                IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN
+           
+                   month_of_year=i_mon
+
+                   day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1))))) 
+
+                   GOTO 38
+
+                ENDIF
+
+             38 ENDDO
+          ENDIF
+       !> --------------------------------------------------------------------------------
+       !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step
+       !>TBD: something to calculate the current year is missing
+       ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
+ 
+          !> calculate month_of_year. TBD: test the condition when day_of_year==31
+ 
+          IF (day_of_year==1) THEN  !> this allows to turn from december to January when passing from a year to another
+ 
+             month_of_year = 1
+        
+          ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN
+
+             month_of_year = month_of_year + 1
+
+          ENDIF
+
+          !> calculate day_of_month
+          IF ( month_of_year == 1 ) THEN
+           
+            day_of_month=day_of_year
+
+          ELSE
+
+            day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1)))))
+
+          ENDIF
+
+
+       ELSE
+
+          !> Condition when date_init is provided but it is given in the wrong format
+          message_string = 'date_init not provided in the namelist or'            //          &
+                              ' given in the wrong format: MUST BE DDMMYYYY'                  
+          CALL message( 'init_date_and_time', 'DT0102', 2, 2, 0, 6, 0 ) 
+
+       ENDIF
+
+
+    END SUBROUTINE init_date_and_time
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate current date and time of the simulation 
+!------------------------------------------------------------------------------!
+ 
+    SUBROUTINE calc_date_and_time
+
+       IMPLICIT NONE
+
+!--    Variables Definition
+       INTEGER                          :: i_mon       !< Index for going through the different months
+
+       !> Update simulation time in seconds
+       time_update = simulated_time-coupling_start_time
+
+!--    Calculate current day of the simulated time
+       days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
+                               / 86400.0_wp ) )
+
+!--    Calculate actual UTC time                        
+       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
+       
+!sB    PRILIMINARY workaround for time_utc changes due to changes in time_since_reference_point in
+!sB    radiation_model_mod during runtime:
+       time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
+
+!--    Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist. 
+
+       IF ( day_of_year_init == 0 ) THEN
+
+          !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM. 
+          IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
+        
+             IF ( month_of_year == 1 ) THEN  !!month of year is read in input
+
+                day_of_year_init = day_of_month
+
+             ELSE
+
+                day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month  !day_of_month is read in input in this case
+
+             ENDIF
+!kanani: Revise, we cannot force users to provide date_init, maybe set a default value?
+!           ELSE
+! 
+!              message_string = 'date_init not provided in the namelist or'             //          &
+!                               ' given in the wrong format: MUST BE DDMMYYYY'                        
+!              CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 ) 
+     
+          ENDIF
+
+       ENDIF
+
       !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59. 
 
@@ -302 +422 @@
 
       ENDIF       
-      
+
     END SUBROUTINE calc_date_and_time
 
-!TBD: check the routines used for update of emission indices in the DEFAULT MODE 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine determines the time factor index in the PRE-PROCESSED emissions mode.
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE time_preprocessed_indices(index_hh)
+
+    USE indices
+
+    IMPLICIT NONE
+
+!
+!-- In/output
+    INTEGER, INTENT(INOUT) ::  index_hh    !> Index Hour
+!
+!-- Additional Variables for calculateing indices
+!-- Constants
+    INTEGER, PARAMETER ::  nhour = 24
+
+    IF (days_since_reference_point == 0) THEN
+
+       index_hh=hour_of_day
+
+    ELSE
+
+       index_hh=(days_since_reference_point*nhour)+(hour_of_day)
+
+    ENDIF
+
+
+ END SUBROUTINE time_preprocessed_indices
+
+
 !------------------------------------------------------------------------------!
 ! Description:

palm/trunk/SOURCE/flow_statistics.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443, basit:
+! bug fixed in chemistry profiles
+! 
+! 3452 2018-10-30 13:13:34Z schwenkel
 ! Bugfix for profiles output
 ! 
@@ -1843 +1847 @@
        ENDIF
 
-       IF ( air_chemistry .AND. max_pr_cs > 0 )  THEN
+       IF ( air_chemistry  .AND.  max_pr_cs > 0 )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-             CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+1,0), sums(nzb,pr_palm+1), &
-                                 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
+             DO  i = 1, max_pr_cs
+                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),       &
+                                    sums(nzb,pr_palm+max_pr_user+i),           &
+                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
+             ENDDO
        ENDIF
 

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443, basit:
+! bug fixed in sums and sums_l for chemistry profile output
+! 
+! 3448 2018-10-29 18:14:31Z kanani
 ! Add biometeorology
 ! 
@@ -711 +715 @@
               ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
               rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
-              sums(nzb:nzt+1,pr_palm+max_pr_user),                             &
-              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),      &
+              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs),                   &
+              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs,0:threads_per_task-1),      &
               sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
               sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                   &

palm/trunk/SOURCE/netcdf_data_input_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443, banzhafs, Russo:
+! Uncommented lines on dimension of surface_fractions
+! Removed par_emis_time_factor variable, moved to chem_emissions_mod
+! Initialized nspec and other emission variables at time of declaration
+! Modified EXPERT mode to PRE-PROCESSED mode
+! Introduced Chemistry static netcdf file
+! Added the routine for reading-in netcdf data for chemistry
+! Added routines to get_variable interface specific for chemistry files
+! 
+! 3429 2018-10-25 13:04:23Z knoop
 ! add default values of origin_x/y/z
 ! 
@@ -368 +378 @@
 
        !-DIMENSIONS
-       INTEGER(iwp)                                 :: nspec                     !< number of chem species for which emission values are provided
-       INTEGER(iwp)                                 :: ncat                      !< number of emission categories
-       INTEGER(iwp)                                 :: nvoc                      !< number of VOCs components
-       INTEGER(iwp)                                 :: npm                       !< number of PMs components
+       INTEGER(iwp)                                 :: nspec=0                   !< number of chem species for which emission values are provided
+       INTEGER(iwp)                                 :: ncat=0                    !< number of emission categories
+       INTEGER(iwp)                                 :: nvoc=0                    !< number of VOCs components
+       INTEGER(iwp)                                 :: npm=0                     !< number of PMs components
        INTEGER(iwp)                                 :: nnox=2                    !< number of NOx components: NO and NO2
        INTEGER(iwp)                                 :: nsox=2                    !< number of SOx components: SO and SO4
@@ -391 +401 @@
        INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: species_index             !< Index of emission chem species
 
-       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),ALLOCATABLE, DIMENSION(:)           :: xm                        !< Molecular masses of emission chem species
 
@@ -928 +935 @@
          ENDDO
 
-          !Assign Constant Values of time factors:
-          emt_att%par_emis_time_factor( : ) = (/ 0.01, 0.01, 0.01, 0.02, 0.03, 0.07, 0.09, 0.09, 0.05, 0.03, 0.03, 0.03, &
-                                                 0.03, 0.03, 0.03, 0.04, 0.05, 0.09, 0.09, 0.09, 0.04, 0.02, 0.01, 0.01 /)
 
        !> DEFAULT AND PRE-PROCESSED MODE

palm/trunk/SOURCE/palm.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! from chemistry branch r3443, forkel:
+! removed double do_emis check around CALL chem_init
+! replaced call to calc_date_and_time to init_date_and_time
+! 
+! 3337 2018-10-12 15:17:09Z kanani
 ! (from branch resler)
 ! Fix chemistry call
@@ -436 +441 @@
 
        CALL chem_init
+
 !       CALL photolysis_init   ! probably also required for restart
-
-       CALL init_date_and_time     !initialize the time of chemistry emissions
 
     ENDIF
@@ -448 +452 @@
 !
 !-- Initialize all necessary variables
-    CALL calc_date_and_time !this is required for chemistry emissions
+!
+!-- Initial time for chem_emissions_mod
+    CALL init_date_and_time
 
     CALL init_3d_model

palm/trunk/SOURCE/prognostic_equations.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! remove duplicate USE chem_modules
+! from chemistry branch r3443, banzhafs, basit:
+! chem_depo call introduced
+! code added for decycling chemistry
+!
+! 3386 2018-10-19 16:28:22Z gronemeier
 ! Renamed tcm_prognostic to tcm_prognostic_equations
 ! 
@@ -340 +346 @@
     USE buoyancy_mod,                                                          &
         ONLY:  buoyancy
-
-    USE chemistry_model_mod,                                                   &
-        ONLY:  chem_integrate, chem_species,  chem_prognostic_equations,       &
-               nspec, nvar, spc_names, chem_boundary_conds
            
     USE chem_modules,                                                          &
-        ONLY:  call_chem_at_all_substeps, chem_gasphase_on
+        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name, do_depo
 
     USE chem_photolysis_mod,                                                   &
         ONLY:  photolysis_control
 
-    USE chem_modules,                                                          &
-        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name
+    USE chemistry_model_mod,                                                   &
+        ONLY:  chem_boundary_conds, chem_depo, chem_integrate,                 &
+               chem_prognostic_equations, chem_species,                        &
+               nspec, nvar, spc_names
 
     USE control_parameters,                                                    &
@@ -495 +499 @@
 
                 IF ( intermediate_timestep_count == 1 .OR.                        &
-                     call_chem_at_all_substeps ) THEN
-                   CALL chem_integrate (i,j)                                                
+                     call_chem_at_all_substeps )  THEN
+                   CALL chem_integrate (i,j)
+                   IF ( do_depo )  THEN
+                      CALL chem_depo(i,j)
+                   ENDIF
                 ENDIF
              ENDDO

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/Makefile

@@ -28 +28 @@
 
 distclean: clean
-        rm -f $(bindir)/kp4.exe
+        rm -f $(bindir)/kpp4palm.exe
 
 # ----------------------------------------------

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh

@@ -40 +40 @@
 # -----------------
 # $Id$ 
+# forkel   30. Oktober 2018: Integrating contents of kp4_compress into chem_gasphase_mod.f90
 # forkel   25. September 2018: Added cat for $MECH to pass mechanism name to kpp4palm
 # ketelsen 18. September 2018: Added cat for '#INLINE F90_GLOBAL' 
@@ -121 +122 @@
       k)   KEEP="YES";;              # keep Working directory
 
-      o)   OUTDIR=$OPTARG;;          # Output directory of Generated Code
-
-      p)   PREFIX=$OPTARG;;          # Name Prefix
-
-      s)   KPP_SOLVER=$OPTARG;;      # Update sample f90 code in the def_MECH directory
-
-      u)   UPDT="YES";;              # keep Working directory
+      o)   OUTDIR=$OPTARG;;          # Output directory of Generated Code = 
+                                     # SOURCE directory (do not change)
+
+      p)   PREFIX=$OPTARG;;          # Name Prefix (chem_gasphase_mod, do not change)
+
+      s)   KPP_SOLVER=$OPTARG;;      # Chosen solver (only Rosebrock solvers work for vector mode)
+
+      u)   UPDT="YES";;              # update mechanisms/def_$MECH/chem_gasphase_mod.f90
 
       v)   MODE="vector";;           # Set to vector Mode
@@ -163 +165 @@
 
 KPP_FILE_LIST="Initialize Integrator LinearAlgebra Jacobian Function Rates Util"
-
+if [[ $MODE = "vector" ]]
+ then
+KPP_FILE_LIST="$KPP_FILE_LIST kp4_compress_subroutines"
+fi
 
 KPP_SUBROUTINE_LIST="Initialize"
@@ -189 +194 @@
    # get vector Solver 
    cp $BASE/templates/${KPP_SOLVER}_vec.f90 ${KPP_SOLVER}.f90
+   cp $BASE/templates/kp4_compress_header ${PREFIX}_kp4_compress_header.f90
+   cp $BASE/templates/kp4_compress_subroutines ${PREFIX}_kp4_compress_subroutines.f90
 fi
 
@@ -207 +214 @@
          KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST Rosenbrock  FunTemplate JacTemplate"
         KPP_INTERFACE_IGNORE="WAXPY"
+
       else 
         KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST FunTemplate JacTemplate"
+        KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST kco_initialize kco_compress kco_finalize"
       fi;;
 
@@ -246 +255 @@
 
 KPP_INCLUDE_LIST="Parameters Global JacobianSP Monitor"
+if [[ $MODE = "vector" ]]
+ then
+KPP_INCLUDE_LIST="$KPP_INCLUDE_LIST kp4_compress_header"
+fi
 
 #Get definition Files
@@ -309 +322 @@
 done
 
-echo start kp4.exe with arguments
+echo start kpp4palm.exe with arguments
 echo $PREFIX $MODE $VLEN $DE_INDEX $DE_INDEX_FAST
 

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C

@@ -12 +12 @@
 //Current revisions:
 //------------------
-/ 
-/ 
+//
+//
 //Former revisions:
 //-----------------
-//$Id: create_kpp_module.C 2470 2017-09-14 13:56:42Z forkel $
-// forkel 25.09.2018:  added  automatic line with meechanism name (read mech_list)
-// forkel 20.09.2018:  added  vl_glo = size(tempi,1)
-// ketelsen 18.09.2018: Removed creation of fill_ Subroutine and creation of calls thereof
-//
-// forkel 12.09.2018: Fix in order not to loose the values of qvap and fakt
-//
-// forkel 03.09.2018: Bug fix: moved kppi.add_line("    CALL initialize after fakt = fakti(is)
-// forkel June 2018: added qvap and fakt, re-established original uppercase/lowercase 
-//                    Deleted definition of qvap,fakt in create_kpp_integrate again
-//
-// ketelsen July 2018: Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate, 
-//                            some cleanup of comments, various changes in create_kpp_integrate)
-//
-// forkel June 2018: Added qvap and fakt, re-established original uppercase/lowercase 
-//
-//
-/ 2017-09-14 13:56:42Z forkel $
-// Removed preprocessor directive __chem again 
-//
-// 2017-09-14 13:56:42Z forkel $
-//
+//$Id: create_kpp_module.C 3453 2018-10-30 13:21:51Z forkel $
+// exclude kco_compress from handling by global_variables2vector (30.10.2018, forkel)
+//
+// Added  automatic line with mechanism name (read mech_list) (25.09.2018, forkel)
+//
+// Added  vl_glo = size(tempi,1) (20.09.2018, forkel) 
+//
+// Removed creation of fill_ Subroutine and creation of calls thereof (18.09.2018, ketelsen)
+//
+// Fix in order not to loose the values of qvap and fakt (12.09.2018, forkel)
+//
+// Bug fixes: moved kppi.add_line("    CALL initialize after fakt = fakti(is)
+// Deleted definition of qvap,fakt in create_kpp_integrate again (03.09.2018, forkel)
+//
+// Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate, 
+// some cleanup of comments, various changes in create_kpp_integrate) (July 2018, ketelsen)
+//
+//
+// Added qvap and fakt                            (June 2018, forkel)
+// --> Change in module_header: qvap, fakt added  (June 2018, forkel)
+//
+// re-established original uppercase/lowercase     (June 2018, forkel)
+// --> Change in module_header: reset case in  Initialize, Integrate, and 
+//     Update_rconst                               (June 2018, forkel)
+//
+// Removed preprocessor directive __chem again (2017-09-14, forkel)
 // 
-//added phot
-//change of some output to lowercase with uppercase Fortran
-//
-// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
-//
-
+// Added phot                                                 (2017-09-14, forkel)
+// --> Change in module_header: Variables for photolyis added (2017-09-14, forkel)
+//
+// change of some output to lowercase with uppercase Fortran (2017, forkel)
+//
+// Intial version of KP4 adapted to PALM                      (Nov. 2016, ketelsen)
+//
 #include <stdio.h>
-//  mz_rs_20090111+
 // stdlib is necessary to define getenv:
 #include <stdlib.h>
-//  mz_rs_20090111-
 
 #include "create_kpp_module.h"
@@ -131 +134 @@
      }
 
-//   Change global variables to vector
+//   Change global variables to vector (except for kp4_compress, which has already the right form)
      
      for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
-       it->global_variables2vector (global_variable_list);
+       if(it->get_name() != "kco_compress" ) {
+         it->global_variables2vector (global_variable_list);
+       }
      }
 
@@ -162 +167 @@
        }
      }
-// cout << "## after Edit individual subroutines    " <<endl;
-
    }
 

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.h

@@ -1 +1 @@
 #ifndef mecca
+#define mecca 1
+//
+// ############################################################################
+//
+//     create_mz_kpp_module
+//
+//     create scalar code from .f90 sources created by KPP to be used in MESSy
+//
+//     COPYRIGHT Klaus Ketelsen and MPI-CH   April 2007
+//
+// ############################################################################
+//
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id$
 // ketelsen 18.09.2018: Removed create_fill_routine
-
-#include <iostream>
-#include <fstream>
-
-#include <string>
-#include <list>
-#include <vector>
-
-#include "fortran_file.h"
-#include "expand_decomp.h"
-
-// Class to create module which contains code generated by kpp.
-
-/ 
-/ 
-// ketelsen 18.09.2018: Removed create_fill_routine
+//
+// initial version                                  (Nov. 2016, ketelsen)
+//
 
 #include <iostream>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C

@@ -9 +9 @@
 //
 // ############################################################################
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: expand_decomp.C 3327 2018-10-09 19:55:00Z forkel $
+//
+// initial version                                  (Nov. 2016, ketelsen)
+
 
 #include <fstream>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.h

@@ -11 +11 @@
 //
 // ############################################################################
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: expand_decomp.h 3327 2018-10-09 19:55:00Z forkel $
+//
+// initial version                                  (Nov. 2016, ketelsen)
+
 
 #include <iostream>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.C

@@ -13 +13 @@
 //Current revisions:
 //------------------
-/ 
-/ 
+// 
+// 
 //Former revisions:
 //-----------------
-//$Id: fortran_file.C 2470 2017-09-14 13:56:42Z forkel $
-//  ketelsen 18.09.2018: (Rev 3260) added vector switch and creation of dimension statement
-//  forkel Sept.2018:  'fill_' routine applied to qvap and fakt
-//                      removal of unnecessary variables (LOOKAT, monitor etc. from Fortran code)
-//  forkel June 2018:  Moved adaption to PALM conventions to the end of the process
-//                     in order to make future use of vector code adaptations possible
-//  forkel 20.04.2018: Replace Roundoff = WLAMCH('E') by Roundoff = epsilon(one)
+//$Id: 
+//
+// Added vector switch and creation of dimension statement (rev. 3260, 18.09.2018, ketelsen)
 // 
+// removal of unnecessary variables (LOOKAT, monitor etc. from Fortran code) (Sept.2018, forkel)
 // 
+// Replace Roundoff = WLAMCH('E') by Roundoff = epsilon(one)   (20.04.2018, forkel)
 // 
-/ 2017-09-14 13:56:42Z forkel $
+// Moved adaption to PALM conventions to the end of the processing kpp output
+//  in order to make future use of vector code adaptations possible (June 2018, forkel)
+// 
+// changed KPP-generated code to lowercase with uppercase Fortran  expressions
+// added photolysis variables                               (2017-09-14, forkel)
 //
-// 
-//changed KPP-generated code to lowercase with uppercase Fortran  expressions
-//added photolysis variables
-//
-//
-//
-//
-// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
+// Initial version of KP4 adapted to PALM                   (Nov 2016, ketelsen)
 //
 

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.h

@@ -14 +14 @@
 //Current revisions:
 //------------------
-/ 
-/ 
-//Former revisions:
-//-----------------
-// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
-//  forkel Sept.2018:  added edit_Initialize
-//  forkel          :  added global_subtolower
 //
-/ Intial version of KP4
-
+//
+// Former revisions:
+// -----------------------
+// $Id: fortran_file.h 3327 2018-10-09 19:55:00Z forkel $
+//
+// added edit_Initialize                            (Sept.2018, forkel)
+//
+// added global_subtolower                          (forkel)
+// initial version                                  (Nov. 2016, ketelsen)
 
 #include <iostream>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file_vec.C

@@ -9 +9 @@
 //
 // ############################################################################
+//
 //Current revisions:
 //------------------
-/ 
-/ 
+//
+//
 //Former revisions:
-//-----------------
-//$Id: fortran_file.C 2470 2017-09-14 13:56:42Z forkel $
-//  ketelsen 18.09.2018: Line 112: do k=is,ie bydo k=1,vl; line 156 ff: replaced index k by j
-//
-//  forkel Sept.2018:  added edit_Initialize
-//                     changed loop direction for update_rconst in edit_Update_RCONST
-/ 2017-09-14 13:56:42Z forkel $
-//
-//
-// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
-//
-
+//-----------------------
+//$Id: fortran_file_vec.C 3327 2018-10-09 19:55:00Z forkel $
+// Line 112: do k=is,ie bydo k=1,vl; line 156 ff: replaced index k by j (18.09.2018, ketelsen)
+//
+// added edit_Initialize, changed loop direction for update_rconst in edit_Update_RCONST (Sept. 2018, forkel)
+//
+// initial version       (Nov. 2016, ketelsen)
+//
 
 #include <fstream>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/main.C

@@ -9 +9 @@
 //
 // ############################################################################
+//
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: main.C 3327 2018-10-09 19:55:00Z forkel $
+//
+// initial version       (Nov. 2016, ketelsen)
+//
 
 //  mz_rs_20090111+

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.C

@@ -9 +9 @@
 //
 // ############################################################################
+//
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: program_line.C 3327 2018-10-09 19:55:00Z forkel $
+//
+// global_subtolower     (June 2018, forkel)
+//
+// initial version       (Nov. 2016, ketelsen)
+//
+
 
 #include "program_line.h"

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.h

@@ -12 +12 @@
 // ############################################################################
 //
-// forkel June 2018: added    void   global_subtolower(string &line);
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: program_line.h 3327 2018-10-09 19:55:00Z forkel $
+//
+// added    void   global_subtolower(string &line); (June 2018, forkel)
+//
+// initial version                                  (Nov. 2016, ketelsen)
 
 #include <iostream>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/utils.C

@@ -9 +9 @@
 //
 // ############################################################################
+//
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: utils.C 3327 2018-10-09 19:55:00Z forkel $
+//
+// initial version                                  (Nov. 2016, ketelsen)
+
 
 //  mz_rs_20090111+

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/utils.h

@@ -11 +11 @@
 //
 // ############################################################################
+//
+//Current revisions:
+//------------------
+//
+//
+// Former revisions:
+// -----------------------
+// $Id: utils.h 3327 2018-10-09 19:55:00Z forkel $
+//
+// initial version                                  (Nov. 2016, ketelsen)
+// 
 
 #include <string>

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90

@@ -2 +2 @@
 !       and can only be used with KP4.
 !
+!Current revisions:
+!------------------
+!
+!
+! Former revisions:
+! -----------------------
+! $Id$
+! commented USE kp4_compress                       (30.10.2018, forkel)
+!
+! initial version (Rev. 3185)                      (June 2018, ketelsen)
+!
+
 SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
            AbsTol,RelTol,              &
@@ -358 +370 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-  USE kp4_compress,    ONLY: kco_initialize, kco_compress, &
-                             kco_finalize, cell_done
+! USE kp4_compress,    ONLY: kco_initialize, kco_compress, &
+!                            kco_finalize, cell_done
   IMPLICIT NONE
 

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header

@@ -40 +40 @@
 !--------------------------------------------------------------------------------!
 !
-! Current revisions:
-! ------------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
-! forkel June 2018: qvap, fakt added
-! forkel June 2018: reset case in  Initialize, Integrate, Update_rconst
-! 
-! 
-! 2460 2017-09-13 14:47:48Z forkel
 !
-! forkel Sept. 2017: Variables for photolyis added
+! MODULE HEADER TEMPLATE
 !
-!
-! Nov. 2016: Intial version (Klaus Ketelsen)
-!
-!------------------------------------------------------------------------------!
-!
-
+!  Initial version (Nov. 2016, ketelsen), for later modifications of module_header 
+!  see comments in kpp4palm/src/create_kpp_module.C
 
 ! Set kpp Double Precision to PALM Default Precision

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_cbm4/chem_gasphase_mod.f90

@@ -44 +44 @@
 !--------------------------------------------------------------------------------!
 !
-! Current revisions:
-! ------------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
-! forkel June 2018: qvap,fakt added
-! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
-! 
-! 
-! 2460 2017-09-13 14:47:48Z forkel
 !
-! forkel Sept. 2017: Variables for photolyis added
+! MODULE HEADER TEMPLATE
 !
-!
-! Nov. 2016: Intial version (Klaus Ketelsen)
-!
-!------------------------------------------------------------------------------!
-!
-
+!  Initial version (Nov. 2016,ketelsen),for later modifications of module_header 
+!  see comments in kpp4palm/src/create_kpp_module.C
 
 ! Set kpp Double Precision to PALM Default Precision
@@ -129 +113 @@
 ! 
 ! File                 : chem_gasphase_mod_Parameters.f90
-! Time                 : Tue Sep 25 18:35:13 2018
-! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Time                 : Tue Oct 30 14:03:29 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
 ! Equation file        : chem_gasphase_mod.kpp
 ! Output root filename : chem_gasphase_mod
@@ -236 +220 @@
 ! 
 ! File                 : chem_gasphase_mod_Global.f90
-! Time                 : Tue Sep 25 18:35:13 2018
-! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Time                 : Tue Oct 30 14:03:29 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20181030mod/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
 ! Equation file        : chem_gasphase_mod.kpp
 ! Output root filename : chem_gasphase_mod