Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3861)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3862)
@@ -27,5 +27,8 @@
 ! -----------------
 ! $Id: chemistry_model_mod.f90 3784 2019-03-05 14:16:20Z banzhafs
-! some formatting 
+! some formatting of the deposition code
+!
+! 3784 2019-03-05 14:16:20Z banzhafs
+! some formatting
 ! 
 ! 3784 2019-03-05 14:16:20Z banzhafs
@@ -2716,5 +2719,5 @@
     IMPLICIT NONE 
 
-    INTEGER(iwp) ::  lsp  !< 
+    INTEGER(iwp) ::  lsp  !< running index for chem spcs. 
 
     DO  lsp = 1, nspec
@@ -2726,4 +2729,43 @@
 
  END SUBROUTINE chem_wrd_local
+
+
+!!! sB remove blanks again later !!!
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
@@ -2766,5 +2808,4 @@
     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
@@ -2784,50 +2825,50 @@
     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))
+!-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
+    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM 
+    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
+    INTEGER(iwp) ::  ind_luv_s_grass = 3                     !<  assigned to ilu_grass
+    INTEGER(iwp) ::  ind_luv_ev_needle_trees = 4             !<  assigned to ilu_coniferous_forest
+    INTEGER(iwp) ::  ind_luv_de_needle_trees = 5             !<  assigned to ilu_coniferous_forest
+    INTEGER(iwp) ::  ind_luv_ev_broad_trees = 6              !<  assigned to ilu_tropical_forest
+    INTEGER(iwp) ::  ind_luv_de_broad_trees = 7              !<  assigned to ilu_deciduous_forest
+    INTEGER(iwp) ::  ind_luv_t_grass = 8                     !<  assigned to ilu_grass
+    INTEGER(iwp) ::  ind_luv_desert = 9                      !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_luv_tundra = 10                     !<  assigned to ilu_other
+    INTEGER(iwp) ::  ind_luv_irr_crops = 11                  !<  assigned to ilu_arable
+    INTEGER(iwp) ::  ind_luv_semidesert = 12                 !<  assigned to ilu_other
+    INTEGER(iwp) ::  ind_luv_ice = 13                        !<  assigned to ilu_ice
+    INTEGER(iwp) ::  ind_luv_marsh = 14                      !<  assigned to ilu_other
+    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
+    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
+    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
+    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to 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 
+    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM  
+    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
+    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
+    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
+    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
+    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to 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
+    INTEGER(iwp) ::  ind_lup_user = 0                        !<  ERROR as no class given in PALM
+    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_asph = 2                        !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert 
+    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert 
+    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_tartan = 13                     !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !<  assigned to ilu_desert
+    INTEGER(iwp) ::  ind_lup_clay = 15                       !<  assigned to ilu_desert
 !
 !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
@@ -2836,20 +2877,20 @@
     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(wp) ::  dt_chem
-    REAL(wp) ::  dh                
-    REAL(wp) ::  inv_dh
-    REAL(wp) ::  dt_dh
+    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
+
+    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
+
+    REAL(wp) ::  dt_chem                !< length of chem time step
+    REAL(wp) ::  dh                     !< vertical grid size 
+    REAL(wp) ::  inv_dh                 !< inverse of vertical grid size 
+    REAL(wp) ::  dt_dh                  !< dt_chem/dh
 
     REAL(wp) ::  dens              !< density at layer k at i,j  
-    REAL(wp) ::  r_aero_lu         !< aerodynamic resistance (s/m) at current surface element
-    REAL(wp) ::  ustar_lu          !< ustar at current surface element
-    REAL(wp) ::  z0h_lu            !< roughness length for heat at current surface element 
-    REAL(wp) ::  glrad             !< rad_sw_in at current surface element 
-    REAL(wp) ::  ppm_to_ugm3       !< conversion factor
-    REAL(wp) ::  relh              !< relative humidity at current surface element
+    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
+    REAL(wp) ::  ustar_surf        !< ustar at current surface element
+    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element 
+    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element 
+    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
+    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
     REAL(wp) ::  lai               !< leaf area index at current surface element
     REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
@@ -2859,24 +2900,24 @@
     REAL(wp) ::  vs                !< Sedimentation velocity
     REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
-    REAL(wp) ::  Rs                !< Sedimentaion resistance (s/m)
-    REAL(wp) ::  Rb                !< quasi-laminar boundary layer resistance (s/m)
-    REAL(wp) ::  Rc_tot            !< total canopy resistance Rc (s/m)
-
-    REAL(wp) ::  catm              !< gasconc. at i, j, k in ug/m3
-    REAL(wp) ::  diffc             !< diffusivity
-
-
-    REAL(wp), DIMENSION(nspec) ::  bud_now_lu       !< budget for LSM vegetation type at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_now_lup      !< budget for LSM pavement type at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_now_luw      !< budget for LSM water type at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_now_luu      !< budget for USM walls/roofs at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_now_lug      !< budget for USM green surfaces at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_now_lud      !< budget for USM windows at current surface element
-    REAL(wp), DIMENSION(nspec) ::  bud_now          !< overall budget at current surface element
-    REAL(wp), DIMENSION(nspec) ::  cc               !< concentration i,j,k
-    REAL(wp), DIMENSION(nspec) ::  ccomp_tot        !< total compensation point (ug/m3)
-    !< For now kept to zero for all species!
-
-    REAL(wp) ::  ttemp          !< temperatur at i,j,k
+    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
+    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
+    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
+
+    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
+    REAL(wp) ::  diffusivity       !< diffusivity
+
+
+    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
+    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
+    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
+    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
+                                                  
+
+    REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
     REAL(wp) ::  ts             !< surface temperatur in degrees celsius
     REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
@@ -3055,8 +3096,8 @@
 !
 !--    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)
+       ustar_surf  = surf_lsm_h%us(m)
+       z0h_surf    = surf_lsm_h%z0h(m)
+       r_aero_surf = surf_lsm_h%r_a(m)
+       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
        lai = surf_lsm_h%lai(m)
        sai = lai + 1
@@ -3064,5 +3105,5 @@
 !--    For small grid spacing neglect R_a
        IF ( dzw(k) <= 1.0 )  THEN
-          r_aero_lu = 0.0_wp
+          r_aero_surf = 0.0_wp
        ENDIF
 !
@@ -3076,49 +3117,49 @@
 ! 
 !--    Initialize budgets
-       bud_now_lu  = 0.0_wp
-       bud_now_lup = 0.0_wp
-       bud_now_luw = 0.0_wp
+       bud_luv = 0.0_wp
+       bud_lup = 0.0_wp
+       bud_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
+          IF ( lu_palm == ind_luv_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
+          ELSEIF ( lu_palm == ind_luv_b_soil )  THEN
              lu_dep = 9
-          ELSEIF ( lu_palm == ind_lu_mixed_crops )  THEN
+          ELSEIF ( lu_palm == ind_luv_mixed_crops )  THEN
              lu_dep = 2
-          ELSEIF ( lu_palm == ind_lu_s_grass )  THEN
+          ELSEIF ( lu_palm == ind_luv_s_grass )  THEN
              lu_dep = 1
-          ELSEIF ( lu_palm == ind_lu_ev_needle_trees )  THEN
+          ELSEIF ( lu_palm == ind_luv_ev_needle_trees )  THEN
              lu_dep = 4
-          ELSEIF ( lu_palm == ind_lu_de_needle_trees )  THEN
+          ELSEIF ( lu_palm == ind_luv_de_needle_trees )  THEN
              lu_dep = 4
-          ELSEIF ( lu_palm == ind_lu_ev_broad_trees )  THEN
+          ELSEIF ( lu_palm == ind_luv_ev_broad_trees )  THEN
              lu_dep = 12
-          ELSEIF ( lu_palm == ind_lu_de_broad_trees )  THEN
+          ELSEIF ( lu_palm == ind_luv_de_broad_trees )  THEN
              lu_dep = 5
-          ELSEIF ( lu_palm == ind_lu_t_grass )  THEN
+          ELSEIF ( lu_palm == ind_luv_t_grass )  THEN
              lu_dep = 1
-          ELSEIF ( lu_palm == ind_lu_desert )  THEN
+          ELSEIF ( lu_palm == ind_luv_desert )  THEN
              lu_dep = 9
-          ELSEIF ( lu_palm == ind_lu_tundra )  THEN
+          ELSEIF ( lu_palm == ind_luv_tundra )  THEN
              lu_dep = 8
-          ELSEIF ( lu_palm == ind_lu_irr_crops )  THEN
+          ELSEIF ( lu_palm == ind_luv_irr_crops )  THEN
              lu_dep = 2
-          ELSEIF ( lu_palm == ind_lu_semidesert )  THEN
+          ELSEIF ( lu_palm == ind_luv_semidesert )  THEN
              lu_dep = 8
-          ELSEIF ( lu_palm == ind_lu_ice )  THEN
+          ELSEIF ( lu_palm == ind_luv_ice )  THEN
              lu_dep = 10
-          ELSEIF ( lu_palm == ind_lu_marsh )  THEN
+          ELSEIF ( lu_palm == ind_luv_marsh )  THEN
              lu_dep = 8
-          ELSEIF ( lu_palm == ind_lu_ev_shrubs )  THEN
+          ELSEIF ( lu_palm == ind_luv_ev_shrubs )  THEN
              lu_dep = 14
-          ELSEIF ( lu_palm == ind_lu_de_shrubs )  THEN
+          ELSEIF ( lu_palm == ind_luv_de_shrubs )  THEN
              lu_dep = 14
-          ELSEIF ( lu_palm == ind_lu_mixed_forest )  THEN
+          ELSEIF ( lu_palm == ind_luv_mixed_forest )  THEN
              lu_dep = 4
-          ELSEIF ( lu_palm == ind_lu_intrup_forest )  THEN
+          ELSEIF ( lu_palm == ind_luv_intrup_forest )  THEN
              lu_dep = 8      
           ENDIF
@@ -3201,14 +3242,14 @@
 !--    Concentration at i,j,k
        DO  lsp = 1, nspec
-          cc(lsp) = chem_species(lsp)%conc(k,j,i)
+          conc_ijk(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
-
-       ts       = ttemp - 273.15  !< in degrees celcius
+       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
+
+       ts       = temp_tmp - 273.15  !< in degrees celcius
 !
 !--    Viscosity of air
-       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
+       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
 !
 !--    Air density at k
@@ -3217,5 +3258,5 @@
 !--    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) )
+       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
 !
 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
@@ -3233,7 +3274,7 @@
              vs = 0.0_wp
              vd_lu = 0.0_wp
-             Rs = 0.0_wp
-             Rb = 0.0_wp
-             Rc_tot = 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
@@ -3245,13 +3286,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      lu_dep,  &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lu(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_luv(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3266,13 +3307,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      lu_dep , &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lu(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_luv(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3302,43 +3343,43 @@
 !--             Use density at k:
 
-                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
+                ppm2ugm3 =  (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
+                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
 !--             Use default value 
-                diffc            = 0.11e-4
+                diffusivity            = 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:
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 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, cos_zenith, &
-                     relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
-                     r_aero_lu , Rb )
+                     z0h_surf, ustar_surf, diffusivity, &
+                     rb )
+!
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
+                     rh_surf, lai, sai, nwet, lu_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
+                     r_aero_surf , rb )
 !
 !--             Calculate budget
-                IF ( Rc_tot <= 0.0 )  THEN
-
-                   bud_now_lu(lsp) = 0.0_wp
+                IF ( rc_tot <= 0.0 )  THEN
+
+                   bud_luv(lsp) = 0.0_wp
 
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-
-                   bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+
+                   bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
@@ -3363,7 +3404,7 @@
              vs = 0.0_wp
              vd_lu = 0.0_wp
-             Rs = 0.0_wp
-             Rb = 0.0_wp
-             Rc_tot = 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
@@ -3375,13 +3416,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      lup_dep,  &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lup(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_lup(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3396,13 +3437,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      lup_dep, &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lup(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_lup(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3433,39 +3474,39 @@
 !--             Use density at lowest layer:
 
-                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
+                ppm2ugm3 =  (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
+                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
 !--             Use default value 
-                diffc            = 0.11e-4
+                diffusivity            = 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:
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 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, cos_zenith, relh, lai, sai, nwet, lup_dep, 2,      &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
-                                         r_aero_lu , Rb )
+                     z0h_surf, ustar_surf, diffusivity, rb )
+!
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
+                                         r_aero_surf , rb )
 !
 !--             Calculate budget
-                IF ( Rc_tot <= 0.0 )  THEN
-                   bud_now_lup(lsp) = 0.0_wp
+                IF ( rc_tot <= 0.0 )  THEN
+                   bud_lup(lsp) = 0.0_wp
                 ELSE
-                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-                   bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+                   bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
@@ -3489,7 +3530,7 @@
              vs = 0.0_wp
              vd_lu = 0.0_wp
-             Rs = 0.0_wp
-             Rb = 0.0_wp
-             Rc_tot = 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
@@ -3501,13 +3542,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      luw_dep, &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_luw(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_luw(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3521,13 +3562,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      luw_dep, &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_luw(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_luw(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3558,43 +3599,43 @@
 !--             Use density at lowest layer:
 
-                ppm_to_ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
+                ppm2ugm3 = (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
+                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
 !--             Use default value 
-                diffc            = 0.11e-4
+                diffusivity            = 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:
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 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, cos_zenith, relh, lai, sai, nwet, luw_dep, 2,      &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
-                                          r_aero_lu , Rb )
+                     z0h_surf, ustar_surf, diffusivity, rb )
+
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         & 
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
+                                         r_aero_surf , rb )
 !
 !--             Calculate budget 
-                IF ( Rc_tot <= 0.0 )  THEN
-
-                   bud_now_luw(lsp) = 0.0_wp
+                IF ( rc_tot <= 0.0 )  THEN
+
+                   bud_luw(lsp) = 0.0_wp
 
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-
-                   bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+
+                   bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
@@ -3605,17 +3646,17 @@
 
 
-       bud_now = 0.0_wp
+       bud = 0.0_wp
 !
 !--    Calculate overall 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)
+          bud(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_luv(lsp) + &
+               surf_lsm_h%frac(ind_pav_green,m) * bud_lup(lsp) + &
+               surf_lsm_h%frac(ind_wat_win,m) * bud_luw(lsp)
 !
 !--       Compute new concentration:
-          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
-
-          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, cc(lsp))
+          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
+
+          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
 
        ENDDO
@@ -3632,8 +3673,8 @@
 !
 !--    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)
+       ustar_surf  = surf_usm_h%us(m)
+       z0h_surf    = surf_usm_h%z0h(m)
+       r_aero_surf = surf_usm_h%r_a(m)
+       solar_rad   = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
        lai = surf_usm_h%lai(m)
        sai = lai + 1
@@ -3641,5 +3682,5 @@
 !--    For small grid spacing neglect R_a
        IF ( dzw(k) <= 1.0 )  THEN
-          r_aero_lu = 0.0_wp
+          r_aero_surf = 0.0_wp
        ENDIF
 !
@@ -3653,7 +3694,7 @@
 !
 !--    Initialize budgets
-       bud_now_luu  = 0.0_wp
-       bud_now_lug = 0.0_wp
-       bud_now_lud = 0.0_wp
+       bud_luu  = 0.0_wp
+       bud_lug = 0.0_wp
+       bud_lud = 0.0_wp
 !
 !--    Get land use for i,j and assign to DEPAC lu
@@ -3662,43 +3703,43 @@
 !--       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
+          lug_palm = ind_luv_s_grass
+          IF ( lug_palm == ind_luv_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
+          ELSEIF ( lug_palm == ind_luv_b_soil )  THEN
              lug_dep = 9
-          ELSEIF ( lug_palm == ind_lu_mixed_crops )  THEN
+          ELSEIF ( lug_palm == ind_luv_mixed_crops )  THEN
              lug_dep = 2
-          ELSEIF ( lug_palm == ind_lu_s_grass )  THEN
+          ELSEIF ( lug_palm == ind_luv_s_grass )  THEN
              lug_dep = 1
-          ELSEIF ( lug_palm == ind_lu_ev_needle_trees )  THEN
+          ELSEIF ( lug_palm == ind_luv_ev_needle_trees )  THEN
              lug_dep = 4
-          ELSEIF ( lug_palm == ind_lu_de_needle_trees )  THEN
+          ELSEIF ( lug_palm == ind_luv_de_needle_trees )  THEN
              lug_dep = 4
-          ELSEIF ( lug_palm == ind_lu_ev_broad_trees )  THEN
+          ELSEIF ( lug_palm == ind_luv_ev_broad_trees )  THEN
              lug_dep = 12
-          ELSEIF ( lug_palm == ind_lu_de_broad_trees )  THEN
+          ELSEIF ( lug_palm == ind_luv_de_broad_trees )  THEN
              lug_dep = 5
-          ELSEIF ( lug_palm == ind_lu_t_grass )  THEN
+          ELSEIF ( lug_palm == ind_luv_t_grass )  THEN
              lug_dep = 1
-          ELSEIF ( lug_palm == ind_lu_desert )  THEN
+          ELSEIF ( lug_palm == ind_luv_desert )  THEN
              lug_dep = 9
-          ELSEIF ( lug_palm == ind_lu_tundra  )  THEN
+          ELSEIF ( lug_palm == ind_luv_tundra  )  THEN
              lug_dep = 8
-          ELSEIF ( lug_palm == ind_lu_irr_crops )  THEN
+          ELSEIF ( lug_palm == ind_luv_irr_crops )  THEN
              lug_dep = 2
-          ELSEIF ( lug_palm == ind_lu_semidesert )  THEN
+          ELSEIF ( lug_palm == ind_luv_semidesert )  THEN
              lug_dep = 8
-          ELSEIF ( lug_palm == ind_lu_ice )  THEN
+          ELSEIF ( lug_palm == ind_luv_ice )  THEN
              lug_dep = 10
-          ELSEIF ( lug_palm == ind_lu_marsh )  THEN
+          ELSEIF ( lug_palm == ind_luv_marsh )  THEN
              lug_dep = 8
-          ELSEIF ( lug_palm == ind_lu_ev_shrubs )  THEN
+          ELSEIF ( lug_palm == ind_luv_ev_shrubs )  THEN
              lug_dep = 14
-          ELSEIF ( lug_palm == ind_lu_de_shrubs  )  THEN
+          ELSEIF ( lug_palm == ind_luv_de_shrubs  )  THEN
              lug_dep = 14
-          ELSEIF ( lug_palm == ind_lu_mixed_forest )  THEN
+          ELSEIF ( lug_palm == ind_luv_mixed_forest )  THEN
              lug_dep = 4
-          ELSEIF ( lug_palm == ind_lu_intrup_forest )  THEN
+          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
              lug_dep = 8      
           ENDIF
@@ -3810,14 +3851,14 @@
 !--    Concentration at i,j,k
        DO  lsp = 1, nspec
-          cc(lsp) = chem_species(lsp)%conc(k,j,i)
+          conc_ijk(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
-
-       ts       = ttemp - 273.15  !< in degrees celcius
+       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
+
+       ts       = temp_tmp - 273.15  !< in degrees celcius
 !
 !--    Viscosity of air
-       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
+       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
 !
 !--    Air density at k
@@ -3826,5 +3867,5 @@
 !--    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) )
+       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
 !
 !--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
@@ -3846,7 +3887,7 @@
              vs = 0.0_wp
              vd_lu = 0.0_wp
-             Rs = 0.0_wp
-             Rb = 0.0_wp
-             Rc_tot = 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
@@ -3858,13 +3899,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      luu_dep,  &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_luu(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_luu(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3878,13 +3919,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      luu_dep , &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_luu(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_luu(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -3914,45 +3955,45 @@
 !--             Use density at k:
 
-                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
+                ppm2ugm3 =  (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
+!--             conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
 !--             Use default value 
-                diffc            = 0.11e-4
+                diffusivity            = 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:
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 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, cos_zenith, relh, lai, sai, nwet, luu_dep, 2,       &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
-                                         r_aero_lu, rb )
+                     z0h_surf, ustar_surf, diffusivity, &
+                     rb )
+!
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,     &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
+                                         r_aero_surf, rb )
 !
 !--             Calculate budget
-                IF ( Rc_tot <= 0.0 )  THEN
-
-                   bud_now_luu(lsp) = 0.0_wp
+                IF ( rc_tot <= 0.0 )  THEN
+
+                   bud_luu(lsp) = 0.0_wp
 
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-
-                   bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+
+                   bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
@@ -3973,7 +4014,7 @@
              vs = 0.0_wp
              vd_lu = 0.0_wp
-             Rs = 0.0_wp
-             Rb = 0.0_wp
-             Rc_tot = 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
@@ -3985,13 +4026,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      lug_dep,  &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lug(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_lug(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -4005,13 +4046,13 @@
                      visc)
 
-                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
+                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, &
+                     nwet, temp_tmp, dens, visc, &
                      lug_dep, &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lug(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_lug(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -4041,44 +4082,44 @@
 !--             Use density at k:
 
-                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
+                ppm2ugm3 =  (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
+                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
 !--             Use default value 
-                diffc            = 0.11e-4
+                diffusivity            = 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:
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 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, cos_zenith, relh, lai, sai, nwet, lug_dep, 2,        &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,              &
-                                         r_aero_lu , Rb )
+                     z0h_surf, ustar_surf, diffusivity, &
+                     rb )
+!
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
+                                         r_aero_surf , rb )
 !
 !--             Calculate budget
-                IF ( Rc_tot <= 0.0 )  THEN
-
-                   bud_now_lug(lsp) = 0.0_wp
+                IF ( rc_tot <= 0.0 )  THEN
+
+                   bud_lug(lsp) = 0.0_wp
 
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-
-                   bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+
+                   bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
@@ -4103,7 +4144,7 @@
              vs = 0.0_wp
              vd_lu = 0.0_wp
-             Rs = 0.0_wp
-             Rb = 0.0_wp
-             Rc_tot = 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
@@ -4118,8 +4159,8 @@
                      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) * &
+                     nwet, temp_tmp, dens, visc,              &
+                     lud_dep, r_aero_surf, ustar_surf )
+
+                bud_lud(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -4136,9 +4177,9 @@
                      particle_pars(ind_p_size, part_type), &
                      particle_pars(ind_p_slip, part_type), &
-                     nwet, ttemp, dens, visc, &
+                     nwet, temp_tmp, dens, visc, &
                      lud_dep, &
-                     r_aero_lu, ustar_lu )
-
-                bud_now_lud(lsp) = - cc(lsp) * &
+                     r_aero_surf, ustar_surf )
+
+                bud_lud(lsp) = - conc_ijk(lsp) * &
                      (1.0_wp - exp(-vd_lu * dt_dh )) * dh
 
@@ -4168,43 +4209,43 @@
 !--             Use density at k:
 
-                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
+                ppm2ugm3 =  (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
+                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
 !
 !--             Diffusivity for DEPAC relevant gases
 !--             Use default value 
-                diffc = 0.11e-4
+                diffusivity = 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:
+                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
+                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
+                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
+                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
+                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
+                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
+                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 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, cos_zenith, relh, lai, sai, nwet, lud_dep, 2,      &
-                                         rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
-                                         r_aero_lu , rb )
+                     z0h_surf, ustar_surf, diffusivity, rb )
+!
+!--             Get rc_tot
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
+                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
+                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
+                                         r_aero_surf , rb )
 !
 !--             Calculate budget 
-                IF ( Rc_tot <= 0.0 )  THEN
-
-                   bud_now_lud(lsp) = 0.0_wp
+                IF ( rc_tot <= 0.0 )  THEN
+
+                   bud_lud(lsp) = 0.0_wp
 
                 ELSE
 
-                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
-
-                   bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
+                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
+
+                   bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
                         (1.0_wp - exp(-vd_lu * dt_dh )) * dh
                 ENDIF
@@ -4215,5 +4256,5 @@
 
 
-       bud_now = 0.0_wp
+       bud = 0.0_wp
 !
 !--    Calculate overall budget for surface m and adapt concentration
@@ -4221,12 +4262,12 @@
 
 
-          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)
+          bud(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_luu(lsp) + &
+               surf_usm_h%frac(ind_pav_green,m) * bud_lug(lsp) + &
+               surf_usm_h%frac(ind_wat_win,m) * bud_lud(lsp)
 !
 !--       Compute new concentration
-          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
-
-          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, cc(lsp) )
+          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
+
+          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
 
        ENDDO
@@ -4250,6 +4291,6 @@
 !> van Zanten et al., 2010.
 !------------------------------------------------------------------------------!
- 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,           &
+ SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
+      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
       ra, rb )  
 !
@@ -4291,12 +4332,12 @@
     REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
     REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
-    REAL(wp), INTENT(IN) ::  glrad                   !< global radiation (W/m2)
+    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
     REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
     REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
     REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
     REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
-    REAL(wp), INTENT(IN) ::  diffc                   !< diffusivity
+    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
     REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
-    REAL(wp), INTENT(IN) ::  catm                    !< actual atmospheric concentration (ug/m3)
+    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
     REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
     REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
@@ -4340,5 +4381,5 @@
 !
 !-- Next statement is just to avoid compiler warning about unused variable
-    IF ( day_of_year == 0  .OR.  ( catm + lat + ra + rb ) > 0.0_wp )  CONTINUE
+    IF ( day_of_year == 0  .OR.  ( conc_ijk_ugm3 + lat + ra + rb ) > 0.0_wp )  CONTINUE
 !
 !-- Define component number
@@ -4406,5 +4447,5 @@
 !
 !--    Stomatal conductance:
-       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
+       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
 !
 !--    Effective soil conductance:
@@ -4419,5 +4460,5 @@
 !--          lai_present, sai_present, &
 !--          ccomp_tot, &
-!--          catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
+!--          conc_ijk_ugm3=conc_ijk_ugm3,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 )
@@ -4426,10 +4467,10 @@
 !--        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'
+!--           IF ( (.not. present(conc_ijk_ugm3)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
+!--              stop 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
 !--           END IF
 !
 !--       Compute rc_eff :
-       !      CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
+       !      CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff)
        !   ENDIF
     ENDIF
@@ -4714,5 +4755,5 @@
 !> Subroutine to compute stomatal conductance 
 !------------------------------------------------------------------------------!
- SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p ) 
+ SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p ) 
 
     IMPLICIT NONE
@@ -4727,9 +4768,9 @@
 
     REAL(wp), INTENT(IN) ::  lai                   !< one-sided leaf area index
-    REAL(wp), INTENT(IN) ::  glrad                 !< global radiation (W/m2)
+    REAL(wp), INTENT(IN) ::  solar_rad             !< solar radiation, dirict+diffuse (W/m2)
     REAL(wp), INTENT(IN) ::  sinphi                !< sin of solar elevation angle
     REAL(wp), INTENT(IN) ::  t                     !< temperature (C)
     REAL(wp), INTENT(IN) ::  rh                    !< relative humidity (%)
-    REAL(wp), INTENT(IN) ::  diffc                 !< diffusion coefficient of the gas involved
+    REAL(wp), INTENT(IN) ::  diffusivity           !< diffusion coefficient of the gas involved
 
     REAL(wp), OPTIONAL,INTENT(IN) :: p             !< pressure (Pa)
@@ -4757,8 +4798,8 @@
        IF ( lai_present )  THEN
 
-          IF ( glrad > 0.0_wp )  THEN
+          IF ( solar_rad > 0.0_wp )  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 
+             CALL rc_gstom_emb( lu, solar_rad, t, vpd, lai_present, lai, sinphi, gstom, p )
+             gstom = gstom * diffusivity / dO3       !< Gstom of Emberson is derived for ozone 
           ELSE
              gstom = 0.0_wp
@@ -4783,5 +4824,5 @@
 !> Subroutine to compute stomatal conductance according to Emberson
 !------------------------------------------------------------------------------!
- SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, lai_present, lai, sinp, Gsto, p )
+ SUBROUTINE rc_gstom_emb( lu, solar_rad, T, vpd, lai_present, lai, sinp, Gsto, p )
 !
 !>  History
@@ -4823,5 +4864,5 @@
     LOGICAL, INTENT(IN) ::  lai_present
 
-    REAL(wp), INTENT(IN) ::  glrad              !< global radiation (W/m2)
+    REAL(wp), INTENT(IN) ::  solar_rad          !< solar radiation, dirict+diffuse (W/m2)
     REAL(wp), INTENT(IN) ::  t                  !< temperature (C)
     REAL(wp), INTENT(IN) ::  vpd                !< vapour pressure deficit (kPa)
@@ -4872,9 +4913,9 @@
 !
 !--    direct and diffuse par, Photoactive (=visible) radiation:
-       CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff )
+       CALL par_dir_diff( solar_rad, 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 > 200.0_wp .AND. lai > 2.5_wp )  THEN
+       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  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
@@ -4883,5 +4924,5 @@
 !--    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 > 200.0_wp .AND. lai > 2.5_wp )  THEN
+       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
           parsun = pardir**0.8 * 0.5 / sinphi + parshade   !< Zhang et al., 2001
        END IF
@@ -4945,13 +4986,13 @@
  !>     Meteorological Service of Canada
  !>     Leiming uses solar irradiance. This should be equal to global radiation and
- !>     Willem Asman set it to global radiation
+ !>     Willem Asman set it to global radiation (here defined as solar radiation, dirict+diffuse)
  !>
  !>     @todo Check/connect/replace with radiation_model_mod variables    
  !-------------------------------------------------------------------
- SUBROUTINE par_dir_diff( glrad, sinphi, pres, pres_0, par_dir, par_diff )
+ SUBROUTINE par_dir_diff( solar_rad, sinphi, pres, pres_0, par_dir, par_diff )
 
     IMPLICIT NONE
 
-    REAL(wp), INTENT(IN) ::  glrad           !< global radiation (W m-2)
+    REAL(wp), INTENT(IN) ::  solar_rad       !< solar radiation, dirict+diffuse (W m-2)
     REAL(wp), INTENT(IN) ::  sinphi          !< sine of the solar elevation
     REAL(wp), INTENT(IN) ::  pres            !< actual pressure (to correct for height) (Pa)
@@ -4995,6 +5036,7 @@
     rn = MAX( 0.01_wp, rdm + rdn )
 !
-!-- Compute ratio between input global radiation and total radiation computed here
-    ratio = MIN( 0.89_wp, glrad / ( rv + rn ) )
+!-- Compute ratio between input global radiation (here defined as solar radiation, dirict+diffuse)
+!-- and total radiation computed here
+    ratio = MIN( 0.89_wp, solar_rad / ( rv + rn ) )
 !
 !-- Calculate total visible radiation
@@ -5015,5 +5057,5 @@
  !> rc_get_vpd: get vapour pressure deficit (kPa)
  !-------------------------------------------------------------------
- SUBROUTINE rc_get_vpd( temp, relh, vpd )
+ SUBROUTINE rc_get_vpd( temp, rh, vpd )
 
     IMPLICIT NONE
@@ -5021,5 +5063,5 @@
 !-- Input/output variables:
     REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
-    REAL(wp), INTENT(IN) ::  relh    !< relative humidity (%)
+    REAL(wp), INTENT(IN) ::  rh    !< relative humidity (%)
 
     REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
@@ -5038,5 +5080,5 @@
 !-- 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 )
+    vpd  = esat * ( 1 - rh / 100 )
 
  END SUBROUTINE rc_get_vpd
@@ -5267,16 +5309,16 @@
  !>
  ! -------------------------------------------------------------------------------------------
-! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, catm, ra, rb, rc_tot, rc_eff )
+! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, conc_ijk_ugm3, ra, rb, rc_tot, rc_eff )
 !
 !    IMPLICIT NONE
 !
 !!-- Input/output variables:
-!    REAL(wp), INTENT(IN) ::  ccomp_tot    !< total compensation point (weighed average of separate compensation points) (ug/m3)
-!    REAL(wp), INTENT(IN) ::  catm         !< atmospheric concentration (ug/m3)
-!    REAL(wp), INTENT(IN) ::  ra           !< aerodynamic resistance (s/m)
-!    REAL(wp), INTENT(IN) ::  rb           !< boundary layer resistance (s/m)
-!    REAL(wp), INTENT(IN) ::  rc_tot       !< total canopy resistance (s/m)
-!
-!    REAL(wp), INTENT(OUT) ::  rc_eff      !< effective total canopy resistance (s/m)
+!    REAL(wp), INTENT(IN) ::  ccomp_tot     !< total compensation point (weighed average of separate compensation points) (ug/m3)
+!    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3 !< atmospheric concentration (ug/m3) above Catm
+!    REAL(wp), INTENT(IN) ::  ra            !< aerodynamic resistance (s/m)
+!    REAL(wp), INTENT(IN) ::  rb            !< boundary layer resistance (s/m)
+!    REAL(wp), INTENT(IN) ::  rc_tot        !< total canopy resistance (s/m)
+!
+!    REAL(wp), INTENT(OUT) ::  rc_eff       !< effective total canopy resistance (s/m)
 !
 !    !
@@ -5287,5 +5329,5 @@
 !       rc_eff = rc_tot
 !
-!    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( catm - ccomp_tot ) < 1.e-8 ) )  THEN
+!    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( conc_ijk_ugm3 - ccomp_tot ) < 1.e-8 ) )  THEN
 !       !
 !!--   surface concentration (almost) equal to atmospheric concentration
@@ -5296,5 +5338,5 @@
 !       !
 !!--    compensation point available, calculate effective resistance
-!       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * catm ) / ( catm - ccomp_tot )
+!       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * conc_ijk_ugm3 ) / ( conc_ijk_ugm3 - ccomp_tot )
 !
 !    ELSE
@@ -5352,6 +5394,6 @@
  !> 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 )
+ SUBROUTINE drydepo_aero_zhang_vd( vd, rs, vs1, partsize, slipcor, nwet, tsurf, dens1, viscos1, &
+      luc, ftop_lu, ustar )
 
     IMPLICIT NONE  
@@ -5370,8 +5412,8 @@
     REAL(wp), INTENT(IN) ::  viscos1         !< air viscosity in lowest layer
     REAL(wp), INTENT(IN) ::  ftop_lu         !< atmospheric resistnace Ra
-    REAL(wp), INTENT(IN) ::  ustar_lu        !< friction velocity u*    
+    REAL(wp), INTENT(IN) ::  ustar           !< friction velocity u*    
 
     REAL(wp), INTENT(OUT) ::  vd             !< deposition velocity (m/s)
-    REAL(wp), INTENT(OUT) ::  Rs             !< sedimentaion resistance (s/m)
+    REAL(wp), INTENT(OUT) ::  rs             !< sedimentaion resistance (s/m)
 !
 !-- constants 
@@ -5417,13 +5459,13 @@
 !-- and sticking efficiency R (1 = no rebound)
     IF ( luc == ilu_ice .OR. nwet==9 .OR. luc == ilu_water_sea .OR. luc == ilu_water_inland )  THEN
-       stokes = vs1 * ustar_lu**2 / ( grav * kinvisc )
+       stokes = vs1 * ustar**2 / ( grav * kinvisc )
        Einterc = 0.0_wp
        Reffic = 1.0_wp
     ELSE IF ( luc == ilu_other .OR. luc == ilu_desert )  THEN     !<tundra of desert
-       stokes = vs1 * ustar_lu**2 / ( grav * kinvisc )
+       stokes = vs1 * ustar**2 / ( grav * kinvisc )
        Einterc = 0.0_wp
        Reffic = exp( -Stokes**0.5_wp )
     ELSE
-       stokes = vs1 * ustar_lu / (grav * A_lu(luc) * 1.e-3)
+       stokes = vs1 * ustar / (grav * A_lu(luc) * 1.e-3)
        Einterc = 0.5_wp * ( partsize / (A_lu(luc) * 1e-3 ) )**2
        Reffic = exp( -Stokes**0.5_wp )
@@ -5437,5 +5479,5 @@
 !
 !-- sedimentation resistance:
-    Rs = 1.0_wp / ( epsilon0 * ustar_lu * ( Ebrown + Eimpac + Einterc ) * Reffic )
+    rs = 1.0_wp / ( epsilon0 * ustar * ( Ebrown + Eimpac + Einterc ) * Reffic )
 
 !-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7):
@@ -5447,5 +5489,5 @@
 !-- where: Rs = Rb (in Seinfeld and Pandis, 2006)
 
-    vd = 1.0_wp / ( ftop_lu + Rs + ftop_lu * Rs * vs1) + vs1
+    vd = 1.0_wp / ( ftop_lu + rs + ftop_lu * rs * vs1) + vs1
 
 
@@ -5457,5 +5499,5 @@
  !> Original EMEP formulation by (Simpson et al, 2003) is used 
  !-------------------------------------------------------------------------------------
- SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb )    
+ SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffusivity, rb )    
 
 
@@ -5469,7 +5511,7 @@
     REAL(wp), INTENT(IN) ::  z0h                  !< roughness length for heat
     REAL(wp), INTENT(IN) ::  ustar                !< friction velocity
-    REAL(wp), INTENT(IN) ::  diffc                !< coefficient of diffusivity
-
-    REAL(wp), INTENT(OUT) ::  Rb                  !< boundary layer resistance
+    REAL(wp), INTENT(IN) ::  diffusivity          !< coefficient of diffusivity
+
+    REAL(wp), INTENT(OUT) ::  rb                  !< boundary layer resistance
 !
 !-- const
@@ -5482,10 +5524,10 @@
 !
 !-- Use Simpson et al. (2003) 
-!-- @TODO: Check Rb over water calculation, until then leave commented lines
+!-- @TODO: Check rb over water calculation, until then leave commented lines
 !--  IF ( is_water )  THEN
-!--   org: Rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.01_wp,ustar))
-!--        Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar))
+!--   org: rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffusivity*kappa_stab*MAX(0.01_wp,ustar))
+!--        rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffusivity*kappa_stab*MAX(0.1_wp,ustar))
 !--  ELSE
-    Rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffc )**0.67_wp
+    rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffusivity )**0.67_wp
 !--  END IF
 
