Ignore:
Timestamp:
Apr 4, 2019 1:07:31 PM (3 years ago)
Author:
banzhafs
Message:

some formatting in deposition code of chemistry_model_mod

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r3848 r3862  
    2727! -----------------
    2828! $Id: chemistry_model_mod.f90 3784 2019-03-05 14:16:20Z banzhafs
    29 ! some formatting
     29! some formatting of the deposition code
     30!
     31! 3784 2019-03-05 14:16:20Z banzhafs
     32! some formatting
    3033!
    3134! 3784 2019-03-05 14:16:20Z banzhafs
     
    27162719    IMPLICIT NONE
    27172720
    2718     INTEGER(iwp) ::  lsp  !<
     2721    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
    27192722
    27202723    DO  lsp = 1, nspec
     
    27262729
    27272730 END SUBROUTINE chem_wrd_local
     2731
     2732
     2733!!! sB remove blanks again later !!!
     2734
     2735
     2736
     2737
     2738
     2739
     2740
     2741
     2742
     2743
     2744
     2745
     2746
     2747
     2748
     2749
     2750
     2751
     2752
     2753
     2754
     2755
     2756
     2757
     2758
     2759
     2760
     2761
     2762
     2763
     2764
     2765
     2766
     2767
     2768
     2769
    27282770
    27292771
     
    27662808    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
    27672809    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
    2768 !    INTEGER(iwp) ::  lu                           !< running index for landuse classes
    27692810    INTEGER(iwp) ::  lu_palm                       !< index of PALM LSM vegetation_type at current surface element
    27702811    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
     
    27842825    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
    27852826!
    2786 !-- Vegetation                                              !<Match to DEPAC
    2787     INTEGER(iwp) ::  ind_lu_user = 0                         !< --> ERROR
    2788     INTEGER(iwp) ::  ind_lu_b_soil = 1                       !< --> ilu_desert
    2789     INTEGER(iwp) ::  ind_lu_mixed_crops = 2                  !< --> ilu_arable
    2790     INTEGER(iwp) ::  ind_lu_s_grass = 3                      !< --> ilu_grass
    2791     INTEGER(iwp) ::  ind_lu_ev_needle_trees = 4              !< --> ilu_coniferous_forest
    2792     INTEGER(iwp) ::  ind_lu_de_needle_trees = 5              !< --> ilu_coniferous_forest
    2793     INTEGER(iwp) ::  ind_lu_ev_broad_trees = 6               !< --> ilu_tropical_forest
    2794     INTEGER(iwp) ::  ind_lu_de_broad_trees = 7               !< --> ilu_deciduous_forest
    2795     INTEGER(iwp) ::  ind_lu_t_grass = 8                      !< --> ilu_grass
    2796     INTEGER(iwp) ::  ind_lu_desert = 9                       !< --> ilu_desert
    2797     INTEGER(iwp) ::  ind_lu_tundra = 10                      !< --> ilu_other
    2798     INTEGER(iwp) ::  ind_lu_irr_crops = 11                   !< --> ilu_arable
    2799     INTEGER(iwp) ::  ind_lu_semidesert = 12                  !< --> ilu_other
    2800     INTEGER(iwp) ::  ind_lu_ice = 13                         !< --> ilu_ice
    2801     INTEGER(iwp) ::  ind_lu_marsh = 14                       !< --> ilu_other
    2802     INTEGER(iwp) ::  ind_lu_ev_shrubs = 15                   !< --> ilu_mediterrean_scrub
    2803     INTEGER(iwp) ::  ind_lu_de_shrubs = 16                   !< --> ilu_mediterrean_scrub
    2804     INTEGER(iwp) ::  ind_lu_mixed_forest = 17                !< --> ilu_coniferous_forest (ave(decid+conif))
    2805     INTEGER(iwp) ::  ind_lu_intrup_forest = 18               !< --> ilu_other (ave(other+decid))
     2827!-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
     2828    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
     2829    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
     2830    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
     2831    INTEGER(iwp) ::  ind_luv_s_grass = 3                     !<  assigned to ilu_grass
     2832    INTEGER(iwp) ::  ind_luv_ev_needle_trees = 4             !<  assigned to ilu_coniferous_forest
     2833    INTEGER(iwp) ::  ind_luv_de_needle_trees = 5             !<  assigned to ilu_coniferous_forest
     2834    INTEGER(iwp) ::  ind_luv_ev_broad_trees = 6              !<  assigned to ilu_tropical_forest
     2835    INTEGER(iwp) ::  ind_luv_de_broad_trees = 7              !<  assigned to ilu_deciduous_forest
     2836    INTEGER(iwp) ::  ind_luv_t_grass = 8                     !<  assigned to ilu_grass
     2837    INTEGER(iwp) ::  ind_luv_desert = 9                      !<  assigned to ilu_desert
     2838    INTEGER(iwp) ::  ind_luv_tundra = 10                     !<  assigned to ilu_other
     2839    INTEGER(iwp) ::  ind_luv_irr_crops = 11                  !<  assigned to ilu_arable
     2840    INTEGER(iwp) ::  ind_luv_semidesert = 12                 !<  assigned to ilu_other
     2841    INTEGER(iwp) ::  ind_luv_ice = 13                        !<  assigned to ilu_ice
     2842    INTEGER(iwp) ::  ind_luv_marsh = 14                      !<  assigned to ilu_other
     2843    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
     2844    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
     2845    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
     2846    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
    28062847!
    28072848!-- Water
    2808     INTEGER(iwp) ::  ind_luw_user = 0                        !< --> ERROR
    2809     INTEGER(iwp) ::  ind_luw_lake = 1                        !< --> ilu_water_inland
    2810     INTEGER(iwp) ::  ind_luw_river = 2                       !< --> ilu_water_inland
    2811     INTEGER(iwp) ::  ind_luw_ocean = 3                       !< --> ilu_water_sea
    2812     INTEGER(iwp) ::  ind_luw_pond = 4                        !< --> ilu_water_inland
    2813     INTEGER(iwp) ::  ind_luw_fountain = 5                    !< --> ilu_water_inland
     2849    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM 
     2850    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
     2851    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
     2852    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
     2853    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
     2854    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
    28142855!
    28152856!-- Pavement
    2816     INTEGER(iwp) ::  ind_lup_user = 0                        !< --> ERROR
    2817     INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !< --> ilu_desert
    2818     INTEGER(iwp) ::  ind_lup_asph = 2                        !< --> ilu_desert
    2819     INTEGER(iwp) ::  ind_lup_conc = 3                        !< --> ilu_desert
    2820     INTEGER(iwp) ::  ind_lup_sett = 4                        !< --> ilu_desert
    2821     INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !< --> ilu_desert
    2822     INTEGER(iwp) ::  ind_lup_cobblest = 6                    !< --> ilu_desert
    2823     INTEGER(iwp) ::  ind_lup_metal = 7                       !< --> ilu_desert
    2824     INTEGER(iwp) ::  ind_lup_wood = 8                        !< --> ilu_desert
    2825     INTEGER(iwp) ::  ind_lup_gravel = 9                      !< --> ilu_desert
    2826     INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !< --> ilu_desert
    2827     INTEGER(iwp) ::  ind_lup_pebblest = 11                   !< --> ilu_desert
    2828     INTEGER(iwp) ::  ind_lup_woodchips = 12                  !< --> ilu_desert
    2829     INTEGER(iwp) ::  ind_lup_tartan = 13                     !< --> ilu_desert
    2830     INTEGER(iwp) ::  ind_lup_art_turf = 14                   !< --> ilu_desert
    2831     INTEGER(iwp) ::  ind_lup_clay = 15                       !< --> ilu_desert
     2857    INTEGER(iwp) ::  ind_lup_user = 0                        !<  ERROR as no class given in PALM
     2858    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !<  assigned to ilu_desert
     2859    INTEGER(iwp) ::  ind_lup_asph = 2                        !<  assigned to ilu_desert
     2860    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
     2861    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
     2862    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
     2863    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
     2864    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
     2865    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
     2866    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
     2867    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
     2868    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
     2869    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !<  assigned to ilu_desert
     2870    INTEGER(iwp) ::  ind_lup_tartan = 13                     !<  assigned to ilu_desert
     2871    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !<  assigned to ilu_desert
     2872    INTEGER(iwp) ::  ind_lup_clay = 15                       !<  assigned to ilu_desert
    28322873!
    28332874!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
     
    28362877    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
    28372878
    2838     INTEGER(iwp) ::  part_type
    2839 
    2840     INTEGER(iwp) ::  nwet           !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
    2841 
    2842     REAL(wp) ::  dt_chem
    2843     REAL(wp) ::  dh               
    2844     REAL(wp) ::  inv_dh
    2845     REAL(wp) ::  dt_dh
     2879    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
     2880
     2881    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     2882
     2883    REAL(wp) ::  dt_chem                !< length of chem time step
     2884    REAL(wp) ::  dh                     !< vertical grid size
     2885    REAL(wp) ::  inv_dh                 !< inverse of vertical grid size
     2886    REAL(wp) ::  dt_dh                  !< dt_chem/dh
    28462887
    28472888    REAL(wp) ::  dens              !< density at layer k at i,j 
    2848     REAL(wp) ::  r_aero_lu         !< aerodynamic resistance (s/m) at current surface element
    2849     REAL(wp) ::  ustar_lu          !< ustar at current surface element
    2850     REAL(wp) ::  z0h_lu            !< roughness length for heat at current surface element
    2851     REAL(wp) ::  glrad             !< rad_sw_in at current surface element
    2852     REAL(wp) ::  ppm_to_ugm3       !< conversion factor
    2853     REAL(wp) ::  relh              !< relative humidity at current surface element
     2889    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
     2890    REAL(wp) ::  ustar_surf        !< ustar at current surface element
     2891    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
     2892    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element
     2893    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
     2894    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
    28542895    REAL(wp) ::  lai               !< leaf area index at current surface element
    28552896    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
     
    28592900    REAL(wp) ::  vs                !< Sedimentation velocity
    28602901    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
    2861     REAL(wp) ::  Rs                !< Sedimentaion resistance (s/m)
    2862     REAL(wp) ::  Rb                !< quasi-laminar boundary layer resistance (s/m)
    2863     REAL(wp) ::  Rc_tot            !< total canopy resistance Rc (s/m)
    2864 
    2865     REAL(wp) ::  catm              !< gasconc. at i, j, k in ug/m3
    2866     REAL(wp) ::  diffc             !< diffusivity
    2867 
    2868 
    2869     REAL(wp), DIMENSION(nspec) ::  bud_now_lu       !< budget for LSM vegetation type at current surface element
    2870     REAL(wp), DIMENSION(nspec) ::  bud_now_lup      !< budget for LSM pavement type at current surface element
    2871     REAL(wp), DIMENSION(nspec) ::  bud_now_luw      !< budget for LSM water type at current surface element
    2872     REAL(wp), DIMENSION(nspec) ::  bud_now_luu      !< budget for USM walls/roofs at current surface element
    2873     REAL(wp), DIMENSION(nspec) ::  bud_now_lug      !< budget for USM green surfaces at current surface element
    2874     REAL(wp), DIMENSION(nspec) ::  bud_now_lud      !< budget for USM windows at current surface element
    2875     REAL(wp), DIMENSION(nspec) ::  bud_now          !< overall budget at current surface element
    2876     REAL(wp), DIMENSION(nspec) ::  cc               !< concentration i,j,k
    2877     REAL(wp), DIMENSION(nspec) ::  ccomp_tot        !< total compensation point (ug/m3)
    2878     !< For now kept to zero for all species!
    2879 
    2880     REAL(wp) ::  ttemp          !< temperatur at i,j,k
     2902    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
     2903    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
     2904    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
     2905
     2906    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
     2907    REAL(wp) ::  diffusivity       !< diffusivity
     2908
     2909
     2910    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
     2911    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
     2912    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
     2913    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
     2914    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
     2915    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
     2916    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
     2917    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
     2918    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
     2919                                                 
     2920
     2921    REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
    28812922    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
    28822923    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
     
    30553096!
    30563097!--    Get needed variables for surface element m
    3057        ustar_lu  = surf_lsm_h%us(m)
    3058        z0h_lu    = surf_lsm_h%z0h(m)
    3059        r_aero_lu = surf_lsm_h%r_a(m)
    3060        glrad     = surf_lsm_h%rad_sw_in(m)
     3098       ustar_surf  = surf_lsm_h%us(m)
     3099       z0h_surf    = surf_lsm_h%z0h(m)
     3100       r_aero_surf = surf_lsm_h%r_a(m)
     3101       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
    30613102       lai = surf_lsm_h%lai(m)
    30623103       sai = lai + 1
     
    30643105!--    For small grid spacing neglect R_a
    30653106       IF ( dzw(k) <= 1.0 )  THEN
    3066           r_aero_lu = 0.0_wp
     3107          r_aero_surf = 0.0_wp
    30673108       ENDIF
    30683109!
     
    30763117!
    30773118!--    Initialize budgets
    3078        bud_now_lu = 0.0_wp
    3079        bud_now_lup = 0.0_wp
    3080        bud_now_luw = 0.0_wp
     3119       bud_luv = 0.0_wp
     3120       bud_lup = 0.0_wp
     3121       bud_luw = 0.0_wp
    30813122!
    30823123!--    Get land use for i,j and assign to DEPAC lu
    30833124       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
    30843125          lu_palm = surf_lsm_h%vegetation_type(m)
    3085           IF ( lu_palm == ind_lu_user )  THEN
     3126          IF ( lu_palm == ind_luv_user )  THEN
    30863127             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
    30873128             CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
    3088           ELSEIF ( lu_palm == ind_lu_b_soil )  THEN
     3129          ELSEIF ( lu_palm == ind_luv_b_soil )  THEN
    30893130             lu_dep = 9
    3090           ELSEIF ( lu_palm == ind_lu_mixed_crops )  THEN
     3131          ELSEIF ( lu_palm == ind_luv_mixed_crops )  THEN
    30913132             lu_dep = 2
    3092           ELSEIF ( lu_palm == ind_lu_s_grass )  THEN
     3133          ELSEIF ( lu_palm == ind_luv_s_grass )  THEN
    30933134             lu_dep = 1
    3094           ELSEIF ( lu_palm == ind_lu_ev_needle_trees )  THEN
     3135          ELSEIF ( lu_palm == ind_luv_ev_needle_trees )  THEN
    30953136             lu_dep = 4
    3096           ELSEIF ( lu_palm == ind_lu_de_needle_trees )  THEN
     3137          ELSEIF ( lu_palm == ind_luv_de_needle_trees )  THEN
    30973138             lu_dep = 4
    3098           ELSEIF ( lu_palm == ind_lu_ev_broad_trees )  THEN
     3139          ELSEIF ( lu_palm == ind_luv_ev_broad_trees )  THEN
    30993140             lu_dep = 12
    3100           ELSEIF ( lu_palm == ind_lu_de_broad_trees )  THEN
     3141          ELSEIF ( lu_palm == ind_luv_de_broad_trees )  THEN
    31013142             lu_dep = 5
    3102           ELSEIF ( lu_palm == ind_lu_t_grass )  THEN
     3143          ELSEIF ( lu_palm == ind_luv_t_grass )  THEN
    31033144             lu_dep = 1
    3104           ELSEIF ( lu_palm == ind_lu_desert )  THEN
     3145          ELSEIF ( lu_palm == ind_luv_desert )  THEN
    31053146             lu_dep = 9
    3106           ELSEIF ( lu_palm == ind_lu_tundra )  THEN
     3147          ELSEIF ( lu_palm == ind_luv_tundra )  THEN
    31073148             lu_dep = 8
    3108           ELSEIF ( lu_palm == ind_lu_irr_crops )  THEN
     3149          ELSEIF ( lu_palm == ind_luv_irr_crops )  THEN
    31093150             lu_dep = 2
    3110           ELSEIF ( lu_palm == ind_lu_semidesert )  THEN
     3151          ELSEIF ( lu_palm == ind_luv_semidesert )  THEN
    31113152             lu_dep = 8
    3112           ELSEIF ( lu_palm == ind_lu_ice )  THEN
     3153          ELSEIF ( lu_palm == ind_luv_ice )  THEN
    31133154             lu_dep = 10
    3114           ELSEIF ( lu_palm == ind_lu_marsh )  THEN
     3155          ELSEIF ( lu_palm == ind_luv_marsh )  THEN
    31153156             lu_dep = 8
    3116           ELSEIF ( lu_palm == ind_lu_ev_shrubs )  THEN
     3157          ELSEIF ( lu_palm == ind_luv_ev_shrubs )  THEN
    31173158             lu_dep = 14
    3118           ELSEIF ( lu_palm == ind_lu_de_shrubs )  THEN
     3159          ELSEIF ( lu_palm == ind_luv_de_shrubs )  THEN
    31193160             lu_dep = 14
    3120           ELSEIF ( lu_palm == ind_lu_mixed_forest )  THEN
     3161          ELSEIF ( lu_palm == ind_luv_mixed_forest )  THEN
    31213162             lu_dep = 4
    3122           ELSEIF ( lu_palm == ind_lu_intrup_forest )  THEN
     3163          ELSEIF ( lu_palm == ind_luv_intrup_forest )  THEN
    31233164             lu_dep = 8     
    31243165          ENDIF
     
    32013242!--    Concentration at i,j,k
    32023243       DO  lsp = 1, nspec
    3203           cc(lsp) = chem_species(lsp)%conc(k,j,i)
     3244          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
    32043245       ENDDO
    32053246
    32063247!--    Temperature at i,j,k
    3207        ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
    3208 
    3209        ts       = ttemp - 273.15  !< in degrees celcius
     3248       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
     3249
     3250       ts       = temp_tmp - 273.15  !< in degrees celcius
    32103251!
    32113252!--    Viscosity of air
    3212        visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
     3253       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
    32133254!
    32143255!--    Air density at k
     
    32173258!--    Calculate relative humidity from specific humidity for DEPAC
    32183259       qv_tmp = MAX(q(k,j,i),0.0_wp)
    3219        relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )
     3260       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
    32203261!
    32213262!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     
    32333274             vs = 0.0_wp
    32343275             vd_lu = 0.0_wp
    3235              Rs = 0.0_wp
    3236              Rb = 0.0_wp
    3237              Rc_tot = 0.0_wp
     3276             rs = 0.0_wp
     3277             rb = 0.0_wp
     3278             rc_tot = 0.0_wp
    32383279             IF ( spc_names(lsp) == 'PM10' )  THEN
    32393280                part_type = 1
     
    32453286                     visc)
    32463287
    3247                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3288                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    32483289                     vs, &
    32493290                     particle_pars(ind_p_size, part_type), &
    32503291                     particle_pars(ind_p_slip, part_type), &
    3251                      nwet, ttemp, dens, visc, &
     3292                     nwet, temp_tmp, dens, visc, &
    32523293                     lu_dep,  &
    3253                      r_aero_lu, ustar_lu )
    3254 
    3255                 bud_now_lu(lsp) = - cc(lsp) * &
     3294                     r_aero_surf, ustar_surf )
     3295
     3296                bud_luv(lsp) = - conc_ijk(lsp) * &
    32563297                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    32573298
     
    32663307                     visc)
    32673308
    3268                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3309                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    32693310                     vs, &
    32703311                     particle_pars(ind_p_size, part_type), &
    32713312                     particle_pars(ind_p_slip, part_type), &
    3272                      nwet, ttemp, dens, visc, &
     3313                     nwet, temp_tmp, dens, visc, &
    32733314                     lu_dep , &
    3274                      r_aero_lu, ustar_lu )
    3275 
    3276                 bud_now_lu(lsp) = - cc(lsp) * &
     3315                     r_aero_surf, ustar_surf )
     3316
     3317                bud_luv(lsp) = - conc_ijk(lsp) * &
    32773318                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    32783319
     
    33023343!--             Use density at k:
    33033344
    3304                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3345                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    33053346!
    33063347!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    33073348                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    3308                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3349                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
    33093350!
    33103351!--             Diffusivity for DEPAC relevant gases
    33113352!--             Use default value
    3312                 diffc            = 0.11e-4
     3353                diffusivity            = 0.11e-4
    33133354!
    33143355!--             overwrite with known coefficients of diffusivity from Massman (1998)
    3315                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    3316                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
    3317                 IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
    3318                 IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
    3319                 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
    3320                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    3321                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3322 !
    3323 !--             Get quasi-laminar boundary layer resistance Rb:
     3356                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     3357                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     3358                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     3359                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     3360                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     3361                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     3362                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     3363!
     3364!--             Get quasi-laminar boundary layer resistance rb:
    33243365                CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), &
    3325                      z0h_lu, ustar_lu, diffc, &
    3326                      Rb )
    3327 !
    3328 !--             Get Rc_tot
    3329                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
    3330                      relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
    3331                      r_aero_lu , Rb )
     3366                     z0h_surf, ustar_surf, diffusivity, &
     3367                     rb )
     3368!
     3369!--             Get rc_tot
     3370                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
     3371                     rh_surf, lai, sai, nwet, lu_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
     3372                     r_aero_surf , rb )
    33323373!
    33333374!--             Calculate budget
    3334                 IF ( Rc_tot <= 0.0 )  THEN
    3335 
    3336                    bud_now_lu(lsp) = 0.0_wp
     3375                IF ( rc_tot <= 0.0 )  THEN
     3376
     3377                   bud_luv(lsp) = 0.0_wp
    33373378
    33383379                ELSE
    33393380
    3340                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    3341 
    3342                    bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3381                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     3382
     3383                   bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    33433384                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    33443385                ENDIF
     
    33633404             vs = 0.0_wp
    33643405             vd_lu = 0.0_wp
    3365              Rs = 0.0_wp
    3366              Rb = 0.0_wp
    3367              Rc_tot = 0.0_wp
     3406             rs = 0.0_wp
     3407             rb = 0.0_wp
     3408             rc_tot = 0.0_wp
    33683409             IF ( spc_names(lsp) == 'PM10' )  THEN
    33693410                part_type = 1
     
    33753416                     visc)
    33763417
    3377                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3418                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    33783419                     vs, &
    33793420                     particle_pars(ind_p_size, part_type), &
    33803421                     particle_pars(ind_p_slip, part_type), &
    3381                      nwet, ttemp, dens, visc, &
     3422                     nwet, temp_tmp, dens, visc, &
    33823423                     lup_dep,  &
    3383                      r_aero_lu, ustar_lu )
    3384 
    3385                 bud_now_lup(lsp) = - cc(lsp) * &
     3424                     r_aero_surf, ustar_surf )
     3425
     3426                bud_lup(lsp) = - conc_ijk(lsp) * &
    33863427                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    33873428
     
    33963437                     visc)
    33973438
    3398                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3439                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    33993440                     vs, &
    34003441                     particle_pars(ind_p_size, part_type), &
    34013442                     particle_pars(ind_p_slip, part_type), &
    3402                      nwet, ttemp, dens, visc, &
     3443                     nwet, temp_tmp, dens, visc, &
    34033444                     lup_dep, &
    3404                      r_aero_lu, ustar_lu )
    3405 
    3406                 bud_now_lup(lsp) = - cc(lsp) * &
     3445                     r_aero_surf, ustar_surf )
     3446
     3447                bud_lup(lsp) = - conc_ijk(lsp) * &
    34073448                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    34083449
     
    34333474!--             Use density at lowest layer:
    34343475
    3435                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3476                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    34363477!
    34373478!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    34383479                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    3439                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3480                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
    34403481!
    34413482!--             Diffusivity for DEPAC relevant gases
    34423483!--             Use default value
    3443                 diffc            = 0.11e-4
     3484                diffusivity            = 0.11e-4
    34443485!
    34453486!--             overwrite with known coefficients of diffusivity from Massman (1998)
    3446                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    3447                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
    3448                 IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
    3449                 IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
    3450                 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
    3451                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    3452                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3453 !
    3454 !--             Get quasi-laminar boundary layer resistance Rb:
     3487                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     3488                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     3489                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     3490                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     3491                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     3492                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     3493                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     3494!
     3495!--             Get quasi-laminar boundary layer resistance rb:
    34553496                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
    3456                      z0h_lu, ustar_lu, diffc, Rb )
    3457 !
    3458 !--             Get Rc_tot
    3459                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
    3460                                          glrad, cos_zenith, relh, lai, sai, nwet, lup_dep, 2,      &
    3461                                          rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
    3462                                          r_aero_lu , Rb )
     3497                     z0h_surf, ustar_surf, diffusivity, rb )
     3498!
     3499!--             Get rc_tot
     3500                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
     3501                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
     3502                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
     3503                                         r_aero_surf , rb )
    34633504!
    34643505!--             Calculate budget
    3465                 IF ( Rc_tot <= 0.0 )  THEN
    3466                    bud_now_lup(lsp) = 0.0_wp
     3506                IF ( rc_tot <= 0.0 )  THEN
     3507                   bud_lup(lsp) = 0.0_wp
    34673508                ELSE
    3468                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    3469                    bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3509                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     3510                   bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    34703511                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    34713512                ENDIF
     
    34893530             vs = 0.0_wp
    34903531             vd_lu = 0.0_wp
    3491              Rs = 0.0_wp
    3492              Rb = 0.0_wp
    3493              Rc_tot = 0.0_wp
     3532             rs = 0.0_wp
     3533             rb = 0.0_wp
     3534             rc_tot = 0.0_wp
    34943535             IF ( spc_names(lsp) == 'PM10' )  THEN
    34953536                part_type = 1
     
    35013542                     visc)
    35023543
    3503                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3544                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    35043545                     vs, &
    35053546                     particle_pars(ind_p_size, part_type), &
    35063547                     particle_pars(ind_p_slip, part_type), &
    3507                      nwet, ttemp, dens, visc, &
     3548                     nwet, temp_tmp, dens, visc, &
    35083549                     luw_dep, &
    3509                      r_aero_lu, ustar_lu )
    3510 
    3511                 bud_now_luw(lsp) = - cc(lsp) * &
     3550                     r_aero_surf, ustar_surf )
     3551
     3552                bud_luw(lsp) = - conc_ijk(lsp) * &
    35123553                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    35133554
     
    35213562                     visc)
    35223563
    3523                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3564                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    35243565                     vs, &
    35253566                     particle_pars(ind_p_size, part_type), &
    35263567                     particle_pars(ind_p_slip, part_type), &
    3527                      nwet, ttemp, dens, visc, &
     3568                     nwet, temp_tmp, dens, visc, &
    35283569                     luw_dep, &
    3529                      r_aero_lu, ustar_lu )
    3530 
    3531                 bud_now_luw(lsp) = - cc(lsp) * &
     3570                     r_aero_surf, ustar_surf )
     3571
     3572                bud_luw(lsp) = - conc_ijk(lsp) * &
    35323573                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    35333574
     
    35583599!--             Use density at lowest layer:
    35593600
    3560                 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3601                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    35613602!
    35623603!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    35633604!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    3564                 catm = cc(lsp) * ppm_to_ugm3 *  specmolm(i_pspec)  ! in ug/m3
     3605                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
    35653606!
    35663607!--             Diffusivity for DEPAC relevant gases
    35673608!--             Use default value
    3568                 diffc            = 0.11e-4
     3609                diffusivity            = 0.11e-4
    35693610!
    35703611!--             overwrite with known coefficients of diffusivity from Massman (1998)
    3571                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    3572                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
    3573                 IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
    3574                 IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
    3575                 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
    3576                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    3577                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3578 !
    3579 !--             Get quasi-laminar boundary layer resistance Rb:
     3612                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     3613                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     3614                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     3615                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     3616                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     3617                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     3618                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     3619!
     3620!--             Get quasi-laminar boundary layer resistance rb:
    35803621                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
    3581                      z0h_lu, ustar_lu, diffc, rb )
    3582 
    3583 !--             Get Rc_tot
    3584                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
    3585                                          glrad, cos_zenith, relh, lai, sai, nwet, luw_dep, 2,      &
    3586                                          rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
    3587                                           r_aero_lu , Rb )
     3622                     z0h_surf, ustar_surf, diffusivity, rb )
     3623
     3624!--             Get rc_tot
     3625                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
     3626                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
     3627                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
     3628                                         r_aero_surf , rb )
    35883629!
    35893630!--             Calculate budget
    3590                 IF ( Rc_tot <= 0.0 )  THEN
    3591 
    3592                    bud_now_luw(lsp) = 0.0_wp
     3631                IF ( rc_tot <= 0.0 )  THEN
     3632
     3633                   bud_luw(lsp) = 0.0_wp
    35933634
    35943635                ELSE
    35953636
    3596                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    3597 
    3598                    bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3637                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     3638
     3639                   bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    35993640                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    36003641                ENDIF
     
    36053646
    36063647
    3607        bud_now = 0.0_wp
     3648       bud = 0.0_wp
    36083649!
    36093650!--    Calculate overall budget for surface m and adapt concentration
    36103651       DO  lsp = 1, nspec
    36113652
    3612           bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_now_lu(lsp) + &
    3613                surf_lsm_h%frac(ind_pav_green,m) * bud_now_lup(lsp) + &
    3614                surf_lsm_h%frac(ind_wat_win,m) * bud_now_luw(lsp)
     3653          bud(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_luv(lsp) + &
     3654               surf_lsm_h%frac(ind_pav_green,m) * bud_lup(lsp) + &
     3655               surf_lsm_h%frac(ind_wat_win,m) * bud_luw(lsp)
    36153656!
    36163657!--       Compute new concentration:
    3617           cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
    3618 
    3619           chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, cc(lsp))
     3658          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
     3659
     3660          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
    36203661
    36213662       ENDDO
     
    36323673!
    36333674!--    Get needed variables for surface element m
    3634        ustar_lu  = surf_usm_h%us(m)
    3635        z0h_lu    = surf_usm_h%z0h(m)
    3636        r_aero_lu = surf_usm_h%r_a(m)
    3637        glrad     = surf_usm_h%rad_sw_in(m)
     3675       ustar_surf  = surf_usm_h%us(m)
     3676       z0h_surf    = surf_usm_h%z0h(m)
     3677       r_aero_surf = surf_usm_h%r_a(m)
     3678       solar_rad   = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
    36383679       lai = surf_usm_h%lai(m)
    36393680       sai = lai + 1
     
    36413682!--    For small grid spacing neglect R_a
    36423683       IF ( dzw(k) <= 1.0 )  THEN
    3643           r_aero_lu = 0.0_wp
     3684          r_aero_surf = 0.0_wp
    36443685       ENDIF
    36453686!
     
    36533694!
    36543695!--    Initialize budgets
    3655        bud_now_luu  = 0.0_wp
    3656        bud_now_lug = 0.0_wp
    3657        bud_now_lud = 0.0_wp
     3696       bud_luu  = 0.0_wp
     3697       bud_lug = 0.0_wp
     3698       bud_lud = 0.0_wp
    36583699!
    36593700!--    Get land use for i,j and assign to DEPAC lu
     
    36623703!--       For green urban surfaces (e.g. green roofs
    36633704!--       assume LU short grass
    3664           lug_palm = ind_lu_s_grass
    3665           IF ( lug_palm == ind_lu_user )  THEN
     3705          lug_palm = ind_luv_s_grass
     3706          IF ( lug_palm == ind_luv_user )  THEN
    36663707             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
    36673708             CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
    3668           ELSEIF ( lug_palm == ind_lu_b_soil )  THEN
     3709          ELSEIF ( lug_palm == ind_luv_b_soil )  THEN
    36693710             lug_dep = 9
    3670           ELSEIF ( lug_palm == ind_lu_mixed_crops )  THEN
     3711          ELSEIF ( lug_palm == ind_luv_mixed_crops )  THEN
    36713712             lug_dep = 2
    3672           ELSEIF ( lug_palm == ind_lu_s_grass )  THEN
     3713          ELSEIF ( lug_palm == ind_luv_s_grass )  THEN
    36733714             lug_dep = 1
    3674           ELSEIF ( lug_palm == ind_lu_ev_needle_trees )  THEN
     3715          ELSEIF ( lug_palm == ind_luv_ev_needle_trees )  THEN
    36753716             lug_dep = 4
    3676           ELSEIF ( lug_palm == ind_lu_de_needle_trees )  THEN
     3717          ELSEIF ( lug_palm == ind_luv_de_needle_trees )  THEN
    36773718             lug_dep = 4
    3678           ELSEIF ( lug_palm == ind_lu_ev_broad_trees )  THEN
     3719          ELSEIF ( lug_palm == ind_luv_ev_broad_trees )  THEN
    36793720             lug_dep = 12
    3680           ELSEIF ( lug_palm == ind_lu_de_broad_trees )  THEN
     3721          ELSEIF ( lug_palm == ind_luv_de_broad_trees )  THEN
    36813722             lug_dep = 5
    3682           ELSEIF ( lug_palm == ind_lu_t_grass )  THEN
     3723          ELSEIF ( lug_palm == ind_luv_t_grass )  THEN
    36833724             lug_dep = 1
    3684           ELSEIF ( lug_palm == ind_lu_desert )  THEN
     3725          ELSEIF ( lug_palm == ind_luv_desert )  THEN
    36853726             lug_dep = 9
    3686           ELSEIF ( lug_palm == ind_lu_tundra  )  THEN
     3727          ELSEIF ( lug_palm == ind_luv_tundra  )  THEN
    36873728             lug_dep = 8
    3688           ELSEIF ( lug_palm == ind_lu_irr_crops )  THEN
     3729          ELSEIF ( lug_palm == ind_luv_irr_crops )  THEN
    36893730             lug_dep = 2
    3690           ELSEIF ( lug_palm == ind_lu_semidesert )  THEN
     3731          ELSEIF ( lug_palm == ind_luv_semidesert )  THEN
    36913732             lug_dep = 8
    3692           ELSEIF ( lug_palm == ind_lu_ice )  THEN
     3733          ELSEIF ( lug_palm == ind_luv_ice )  THEN
    36933734             lug_dep = 10
    3694           ELSEIF ( lug_palm == ind_lu_marsh )  THEN
     3735          ELSEIF ( lug_palm == ind_luv_marsh )  THEN
    36953736             lug_dep = 8
    3696           ELSEIF ( lug_palm == ind_lu_ev_shrubs )  THEN
     3737          ELSEIF ( lug_palm == ind_luv_ev_shrubs )  THEN
    36973738             lug_dep = 14
    3698           ELSEIF ( lug_palm == ind_lu_de_shrubs  )  THEN
     3739          ELSEIF ( lug_palm == ind_luv_de_shrubs  )  THEN
    36993740             lug_dep = 14
    3700           ELSEIF ( lug_palm == ind_lu_mixed_forest )  THEN
     3741          ELSEIF ( lug_palm == ind_luv_mixed_forest )  THEN
    37013742             lug_dep = 4
    3702           ELSEIF ( lug_palm == ind_lu_intrup_forest )  THEN
     3743          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
    37033744             lug_dep = 8     
    37043745          ENDIF
     
    38103851!--    Concentration at i,j,k
    38113852       DO  lsp = 1, nspec
    3812           cc(lsp) = chem_species(lsp)%conc(k,j,i)
     3853          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
    38133854       ENDDO
    38143855!
    38153856!--    Temperature at i,j,k
    3816        ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
    3817 
    3818        ts       = ttemp - 273.15  !< in degrees celcius
     3857       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
     3858
     3859       ts       = temp_tmp - 273.15  !< in degrees celcius
    38193860!
    38203861!--    Viscosity of air
    3821        visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
     3862       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
    38223863!
    38233864!--    Air density at k
     
    38263867!--    Calculate relative humidity from specific humidity for DEPAC
    38273868       qv_tmp = MAX(q(k,j,i),0.0_wp)
    3828        relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )
     3869       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
    38293870!
    38303871!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     
    38463887             vs = 0.0_wp
    38473888             vd_lu = 0.0_wp
    3848              Rs = 0.0_wp
    3849              Rb = 0.0_wp
    3850              Rc_tot = 0.0_wp
     3889             rs = 0.0_wp
     3890             rb = 0.0_wp
     3891             rc_tot = 0.0_wp
    38513892             IF (spc_names(lsp) == 'PM10' )  THEN
    38523893                part_type = 1
     
    38583899                     visc)
    38593900
    3860                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3901                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    38613902                     vs, &
    38623903                     particle_pars(ind_p_size, part_type), &
    38633904                     particle_pars(ind_p_slip, part_type), &
    3864                      nwet, ttemp, dens, visc, &
     3905                     nwet, temp_tmp, dens, visc, &
    38653906                     luu_dep,  &
    3866                      r_aero_lu, ustar_lu )
    3867 
    3868                 bud_now_luu(lsp) = - cc(lsp) * &
     3907                     r_aero_surf, ustar_surf )
     3908
     3909                bud_luu(lsp) = - conc_ijk(lsp) * &
    38693910                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    38703911
     
    38783919                     visc)
    38793920
    3880                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3921                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    38813922                     vs, &
    38823923                     particle_pars(ind_p_size, part_type), &
    38833924                     particle_pars(ind_p_slip, part_type), &
    3884                      nwet, ttemp, dens, visc, &
     3925                     nwet, temp_tmp, dens, visc, &
    38853926                     luu_dep , &
    3886                      r_aero_lu, ustar_lu )
    3887 
    3888                 bud_now_luu(lsp) = - cc(lsp) * &
     3927                     r_aero_surf, ustar_surf )
     3928
     3929                bud_luu(lsp) = - conc_ijk(lsp) * &
    38893930                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    38903931
     
    39143955!--             Use density at k:
    39153956
    3916                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3957                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
    39173958
    39183959                !
    39193960!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    39203961!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    3921 !--             catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3962!--             conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
    39223963!
    39233964!--             Diffusivity for DEPAC relevant gases
    39243965!--             Use default value
    3925                 diffc            = 0.11e-4
     3966                diffusivity            = 0.11e-4
    39263967!
    39273968!--             overwrite with known coefficients of diffusivity from Massman (1998)
    3928                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    3929                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
    3930                 IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
    3931                 IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
    3932                 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
    3933                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    3934                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    3935 !
    3936 !--             Get quasi-laminar boundary layer resistance Rb:
     3969                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     3970                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     3971                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     3972                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     3973                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     3974                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     3975                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     3976!
     3977!--             Get quasi-laminar boundary layer resistance rb:
    39373978                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
    3938                      z0h_lu, ustar_lu, diffc, &
    3939                      Rb )
    3940 !
    3941 !--             Get Rc_tot
    3942                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,       &
    3943                                          glrad, cos_zenith, relh, lai, sai, nwet, luu_dep, 2,       &
    3944                                          rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,             &
    3945                                          r_aero_lu, rb )
     3979                     z0h_surf, ustar_surf, diffusivity, &
     3980                     rb )
     3981!
     3982!--             Get rc_tot
     3983                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
     3984                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,     &
     3985                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
     3986                                         r_aero_surf, rb )
    39463987!
    39473988!--             Calculate budget
    3948                 IF ( Rc_tot <= 0.0 )  THEN
    3949 
    3950                    bud_now_luu(lsp) = 0.0_wp
     3989                IF ( rc_tot <= 0.0 )  THEN
     3990
     3991                   bud_luu(lsp) = 0.0_wp
    39513992
    39523993                ELSE
    39533994
    3954                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    3955 
    3956                    bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3995                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     3996
     3997                   bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    39573998                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    39583999                ENDIF
     
    39734014             vs = 0.0_wp
    39744015             vd_lu = 0.0_wp
    3975              Rs = 0.0_wp
    3976              Rb = 0.0_wp
    3977              Rc_tot = 0.0_wp
     4016             rs = 0.0_wp
     4017             rb = 0.0_wp
     4018             rc_tot = 0.0_wp
    39784019             IF ( spc_names(lsp) == 'PM10' )  THEN
    39794020                part_type = 1
     
    39854026                     visc)
    39864027
    3987                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     4028                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    39884029                     vs, &
    39894030                     particle_pars(ind_p_size, part_type), &
    39904031                     particle_pars(ind_p_slip, part_type), &
    3991                      nwet, ttemp, dens, visc, &
     4032                     nwet, temp_tmp, dens, visc, &
    39924033                     lug_dep,  &
    3993                      r_aero_lu, ustar_lu )
    3994 
    3995                 bud_now_lug(lsp) = - cc(lsp) * &
     4034                     r_aero_surf, ustar_surf )
     4035
     4036                bud_lug(lsp) = - conc_ijk(lsp) * &
    39964037                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    39974038
     
    40054046                     visc)
    40064047
    4007                 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     4048                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
    40084049                     vs, &
    40094050                     particle_pars(ind_p_size, part_type), &
    40104051                     particle_pars(ind_p_slip, part_type), &
    4011                      nwet, ttemp, dens, visc, &
     4052                     nwet, temp_tmp, dens, visc, &
    40124053                     lug_dep, &
    4013                      r_aero_lu, ustar_lu )
    4014 
    4015                 bud_now_lug(lsp) = - cc(lsp) * &
     4054                     r_aero_surf, ustar_surf )
     4055
     4056                bud_lug(lsp) = - conc_ijk(lsp) * &
    40164057                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    40174058
     
    40414082!--             Use density at k:
    40424083
    4043                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
     4084                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
    40444085!
    40454086!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    40464087                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    4047                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     4088                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
    40484089!
    40494090!--             Diffusivity for DEPAC relevant gases
    40504091!--             Use default value
    4051                 diffc            = 0.11e-4
     4092                diffusivity            = 0.11e-4
    40524093!
    40534094!--             overwrite with known coefficients of diffusivity from Massman (1998)
    4054                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    4055                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
    4056                 IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
    4057                 IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
    4058                 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
    4059                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    4060                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    4061 !
    4062 !--             Get quasi-laminar boundary layer resistance Rb:
     4095                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     4096                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     4097                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     4098                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     4099                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     4100                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     4101                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     4102!
     4103!--             Get quasi-laminar boundary layer resistance rb:
    40634104                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
    4064                      z0h_lu, ustar_lu, diffc, &
    4065                      Rb )
    4066 !
    4067 !--             Get Rc_tot
    4068                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,        &
    4069                                          glrad, cos_zenith, relh, lai, sai, nwet, lug_dep, 2,        &
    4070                                          rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,              &
    4071                                          r_aero_lu , Rb )
     4105                     z0h_surf, ustar_surf, diffusivity, &
     4106                     rb )
     4107!
     4108!--             Get rc_tot
     4109                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
     4110                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
     4111                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
     4112                                         r_aero_surf , rb )
    40724113!
    40734114!--             Calculate budget
    4074                 IF ( Rc_tot <= 0.0 )  THEN
    4075 
    4076                    bud_now_lug(lsp) = 0.0_wp
     4115                IF ( rc_tot <= 0.0 )  THEN
     4116
     4117                   bud_lug(lsp) = 0.0_wp
    40774118
    40784119                ELSE
    40794120
    4080                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    4081 
    4082                    bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     4121                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     4122
     4123                   bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    40834124                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    40844125                ENDIF
     
    41034144             vs = 0.0_wp
    41044145             vd_lu = 0.0_wp
    4105              Rs = 0.0_wp
    4106              Rb = 0.0_wp
    4107              Rc_tot = 0.0_wp
     4146             rs = 0.0_wp
     4147             rb = 0.0_wp
     4148             rc_tot = 0.0_wp
    41084149             IF ( spc_names(lsp) == 'PM10' )  THEN
    41094150                part_type = 1
     
    41184159                     particle_pars(ind_p_size, part_type), &
    41194160                     particle_pars(ind_p_slip, part_type), &
    4120                      nwet, ttemp, dens, visc,              &
    4121                      lud_dep, r_aero_lu, ustar_lu )
    4122 
    4123                 bud_now_lud(lsp) = - cc(lsp) * &
     4161                     nwet, temp_tmp, dens, visc,              &
     4162                     lud_dep, r_aero_surf, ustar_surf )
     4163
     4164                bud_lud(lsp) = - conc_ijk(lsp) * &
    41244165                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    41254166
     
    41364177                     particle_pars(ind_p_size, part_type), &
    41374178                     particle_pars(ind_p_slip, part_type), &
    4138                      nwet, ttemp, dens, visc, &
     4179                     nwet, temp_tmp, dens, visc, &
    41394180                     lud_dep, &
    4140                      r_aero_lu, ustar_lu )
    4141 
    4142                 bud_now_lud(lsp) = - cc(lsp) * &
     4181                     r_aero_surf, ustar_surf )
     4182
     4183                bud_lud(lsp) = - conc_ijk(lsp) * &
    41434184                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    41444185
     
    41684209!--             Use density at k:
    41694210
    4170                 ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
     4211                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
    41714212!
    41724213!--             Atmospheric concentration in DEPAC is requested in ug/m3:
    41734214!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
    4174                 catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     4215                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
    41754216!
    41764217!--             Diffusivity for DEPAC relevant gases
    41774218!--             Use default value
    4178                 diffc = 0.11e-4
     4219                diffusivity = 0.11e-4
    41794220!
    41804221!--             overwrite with known coefficients of diffusivity from Massman (1998)
    4181                 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
    4182                 IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
    4183                 IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
    4184                 IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
    4185                 IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
    4186                 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
    4187                 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
    4188 !
    4189 !--             Get quasi-laminar boundary layer resistance Rb:
     4222                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
     4223                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
     4224                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
     4225                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
     4226                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
     4227                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
     4228                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
     4229!
     4230!--             Get quasi-laminar boundary layer resistance rb:
    41904231                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
    4191                      z0h_lu, ustar_lu, diffc, rb )
    4192 !
    4193 !--             Get Rc_tot
    4194                 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu,      &
    4195                                          glrad, cos_zenith, relh, lai, sai, nwet, lud_dep, 2,      &
    4196                                          rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc,            &
    4197                                          r_aero_lu , rb )
     4232                     z0h_surf, ustar_surf, diffusivity, rb )
     4233!
     4234!--             Get rc_tot
     4235                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
     4236                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
     4237                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
     4238                                         r_aero_surf , rb )
    41984239!
    41994240!--             Calculate budget
    4200                 IF ( Rc_tot <= 0.0 )  THEN
    4201 
    4202                    bud_now_lud(lsp) = 0.0_wp
     4241                IF ( rc_tot <= 0.0 )  THEN
     4242
     4243                   bud_lud(lsp) = 0.0_wp
    42034244
    42044245                ELSE
    42054246
    4206                    vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
    4207 
    4208                    bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     4247                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
     4248
     4249                   bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
    42094250                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
    42104251                ENDIF
     
    42154256
    42164257
    4217        bud_now = 0.0_wp
     4258       bud = 0.0_wp
    42184259!
    42194260!--    Calculate overall budget for surface m and adapt concentration
     
    42214262
    42224263
    4223           bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_now_luu(lsp) + &
    4224                surf_usm_h%frac(ind_pav_green,m) * bud_now_lug(lsp) + &
    4225                surf_usm_h%frac(ind_wat_win,m) * bud_now_lud(lsp)
     4264          bud(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_luu(lsp) + &
     4265               surf_usm_h%frac(ind_pav_green,m) * bud_lug(lsp) + &
     4266               surf_usm_h%frac(ind_wat_win,m) * bud_lud(lsp)
    42264267!
    42274268!--       Compute new concentration
    4228           cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
    4229 
    4230           chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, cc(lsp) )
     4269          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
     4270
     4271          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
    42314272
    42324273       ENDDO
     
    42504291!> van Zanten et al., 2010.
    42514292!------------------------------------------------------------------------------!
    4252  SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, glrad, sinphi,  &
    4253       rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc,           &
     4293 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
     4294      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
    42544295      ra, rb ) 
    42554296!
     
    42914332    REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
    42924333    REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
    4293     REAL(wp), INTENT(IN) ::  glrad                   !< global radiation (W/m2)
     4334    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
    42944335    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
    42954336    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
    42964337    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
    42974338    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
    4298     REAL(wp), INTENT(IN) ::  diffc                   !< diffusivity
     4339    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
    42994340    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
    4300     REAL(wp), INTENT(IN) ::  catm                    !< actual atmospheric concentration (ug/m3)
     4341    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
    43014342    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
    43024343    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
     
    43404381!
    43414382!-- Next statement is just to avoid compiler warning about unused variable
    4342     IF ( day_of_year == 0  .OR.  ( catm + lat + ra + rb ) > 0.0_wp )  CONTINUE
     4383    IF ( day_of_year == 0  .OR.  ( conc_ijk_ugm3 + lat + ra + rb ) > 0.0_wp )  CONTINUE
    43434384!
    43444385!-- Define component number
     
    44064447!
    44074448!--    Stomatal conductance:
    4408        CALL rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
     4449       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
    44094450!
    44104451!--    Effective soil conductance:
     
    44194460!--          lai_present, sai_present, &
    44204461!--          ccomp_tot, &
    4421 !--          catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
     4462!--          conc_ijk_ugm3=conc_ijk_ugm3,c_ave_prev_nh3=c_ave_prev_nh3, &
    44224463!--          c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
    44234464!--          tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
     
    44264467!--        IF ( present(rc_eff) ) then
    44274468!--          check on required arguments:
    4428 !--           IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
    4429 !--              stop 'output argument rc_eff requires input arguments catm, ra and rb'
     4469!--           IF ( (.not. present(conc_ijk_ugm3)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
     4470!--              stop 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
    44304471!--           END IF
    44314472!
    44324473!--       Compute rc_eff :
    4433        !      CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
     4474       !      CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff)
    44344475       !   ENDIF
    44354476    ENDIF
     
    47144755!> Subroutine to compute stomatal conductance
    47154756!------------------------------------------------------------------------------!
    4716  SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
     4757 SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
    47174758
    47184759    IMPLICIT NONE
     
    47274768
    47284769    REAL(wp), INTENT(IN) ::  lai                   !< one-sided leaf area index
    4729     REAL(wp), INTENT(IN) ::  glrad                 !< global radiation (W/m2)
     4770    REAL(wp), INTENT(IN) ::  solar_rad             !< solar radiation, dirict+diffuse (W/m2)
    47304771    REAL(wp), INTENT(IN) ::  sinphi                !< sin of solar elevation angle
    47314772    REAL(wp), INTENT(IN) ::  t                     !< temperature (C)
    47324773    REAL(wp), INTENT(IN) ::  rh                    !< relative humidity (%)
    4733     REAL(wp), INTENT(IN) ::  diffc                 !< diffusion coefficient of the gas involved
     4774    REAL(wp), INTENT(IN) ::  diffusivity           !< diffusion coefficient of the gas involved
    47344775
    47354776    REAL(wp), OPTIONAL,INTENT(IN) :: p             !< pressure (Pa)
     
    47574798       IF ( lai_present )  THEN
    47584799
    4759           IF ( glrad > 0.0_wp )  THEN
     4800          IF ( solar_rad > 0.0_wp )  THEN
    47604801             CALL rc_get_vpd( t, rh, vpd )
    4761              CALL rc_gstom_emb( lu, glrad, t, vpd, lai_present, lai, sinphi, gstom, p )
    4762              gstom = gstom * diffc / dO3       !< Gstom of Emberson is derived for ozone
     4802             CALL rc_gstom_emb( lu, solar_rad, t, vpd, lai_present, lai, sinphi, gstom, p )
     4803             gstom = gstom * diffusivity / dO3       !< Gstom of Emberson is derived for ozone
    47634804          ELSE
    47644805             gstom = 0.0_wp
     
    47834824!> Subroutine to compute stomatal conductance according to Emberson
    47844825!------------------------------------------------------------------------------!
    4785  SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, lai_present, lai, sinp, Gsto, p )
     4826 SUBROUTINE rc_gstom_emb( lu, solar_rad, T, vpd, lai_present, lai, sinp, Gsto, p )
    47864827!
    47874828!>  History
     
    48234864    LOGICAL, INTENT(IN) ::  lai_present
    48244865
    4825     REAL(wp), INTENT(IN) ::  glrad              !< global radiation (W/m2)
     4866    REAL(wp), INTENT(IN) ::  solar_rad          !< solar radiation, dirict+diffuse (W/m2)
    48264867    REAL(wp), INTENT(IN) ::  t                  !< temperature (C)
    48274868    REAL(wp), INTENT(IN) ::  vpd                !< vapour pressure deficit (kPa)
     
    48724913!
    48734914!--    direct and diffuse par, Photoactive (=visible) radiation:
    4874        CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff )
     4915       CALL par_dir_diff( solar_rad, sinphi, pres, p_sealevel, pardir, pardiff )
    48754916!
    48764917!--    par for shaded leaves (canopy averaged):
    48774918       parshade = pardiff * exp( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )     !< Norman,1982
    4878        IF ( glrad > 200.0_wp .AND. lai > 2.5_wp )  THEN
     4919       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
    48794920          parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )  !< Zhang et al., 2001
    48804921       END IF
     
    48834924!--    alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
    48844925       parsun = pardir * 0.5/sinphi + parshade             !< Norman, 1982
    4885        IF ( glrad > 200.0_wp .AND. lai > 2.5_wp )  THEN
     4926       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
    48864927          parsun = pardir**0.8 * 0.5 / sinphi + parshade   !< Zhang et al., 2001
    48874928       END IF
     
    49454986 !>     Meteorological Service of Canada
    49464987 !>     Leiming uses solar irradiance. This should be equal to global radiation and
    4947  !>     Willem Asman set it to global radiation
     4988 !>     Willem Asman set it to global radiation (here defined as solar radiation, dirict+diffuse)
    49484989 !>
    49494990 !>     @todo Check/connect/replace with radiation_model_mod variables   
    49504991 !-------------------------------------------------------------------
    4951  SUBROUTINE par_dir_diff( glrad, sinphi, pres, pres_0, par_dir, par_diff )
     4992 SUBROUTINE par_dir_diff( solar_rad, sinphi, pres, pres_0, par_dir, par_diff )
    49524993
    49534994    IMPLICIT NONE
    49544995
    4955     REAL(wp), INTENT(IN) ::  glrad           !< global radiation (W m-2)
     4996    REAL(wp), INTENT(IN) ::  solar_rad       !< solar radiation, dirict+diffuse (W m-2)
    49564997    REAL(wp), INTENT(IN) ::  sinphi          !< sine of the solar elevation
    49574998    REAL(wp), INTENT(IN) ::  pres            !< actual pressure (to correct for height) (Pa)
     
    49955036    rn = MAX( 0.01_wp, rdm + rdn )
    49965037!
    4997 !-- Compute ratio between input global radiation and total radiation computed here
    4998     ratio = MIN( 0.89_wp, glrad / ( rv + rn ) )
     5038!-- Compute ratio between input global radiation (here defined as solar radiation, dirict+diffuse)
     5039!-- and total radiation computed here
     5040    ratio = MIN( 0.89_wp, solar_rad / ( rv + rn ) )
    49995041!
    50005042!-- Calculate total visible radiation
     
    50155057 !> rc_get_vpd: get vapour pressure deficit (kPa)
    50165058 !-------------------------------------------------------------------
    5017  SUBROUTINE rc_get_vpd( temp, relh, vpd )
     5059 SUBROUTINE rc_get_vpd( temp, rh, vpd )
    50185060
    50195061    IMPLICIT NONE
     
    50215063!-- Input/output variables:
    50225064    REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
    5023     REAL(wp), INTENT(IN) ::  relh    !< relative humidity (%)
     5065    REAL(wp), INTENT(IN) ::  rh    !< relative humidity (%)
    50245066
    50255067    REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
     
    50385080!-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973)
    50395081    esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5
    5040     vpd  = esat * ( 1 - relh / 100 )
     5082    vpd  = esat * ( 1 - rh / 100 )
    50415083
    50425084 END SUBROUTINE rc_get_vpd
     
    52675309 !>
    52685310 ! -------------------------------------------------------------------------------------------
    5269 ! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, catm, ra, rb, rc_tot, rc_eff )
     5311! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, conc_ijk_ugm3, ra, rb, rc_tot, rc_eff )
    52705312!
    52715313!    IMPLICIT NONE
    52725314!
    52735315!!-- Input/output variables:
    5274 !    REAL(wp), INTENT(IN) ::  ccomp_tot    !< total compensation point (weighed average of separate compensation points) (ug/m3)
    5275 !    REAL(wp), INTENT(IN) ::  catm         !< atmospheric concentration (ug/m3)
    5276 !    REAL(wp), INTENT(IN) ::  ra           !< aerodynamic resistance (s/m)
    5277 !    REAL(wp), INTENT(IN) ::  rb           !< boundary layer resistance (s/m)
    5278 !    REAL(wp), INTENT(IN) ::  rc_tot       !< total canopy resistance (s/m)
    5279 !
    5280 !    REAL(wp), INTENT(OUT) ::  rc_eff      !< effective total canopy resistance (s/m)
     5316!    REAL(wp), INTENT(IN) ::  ccomp_tot     !< total compensation point (weighed average of separate compensation points) (ug/m3)
     5317!    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3 !< atmospheric concentration (ug/m3) above Catm
     5318!    REAL(wp), INTENT(IN) ::  ra            !< aerodynamic resistance (s/m)
     5319!    REAL(wp), INTENT(IN) ::  rb            !< boundary layer resistance (s/m)
     5320!    REAL(wp), INTENT(IN) ::  rc_tot        !< total canopy resistance (s/m)
     5321!
     5322!    REAL(wp), INTENT(OUT) ::  rc_eff       !< effective total canopy resistance (s/m)
    52815323!
    52825324!    !
     
    52875329!       rc_eff = rc_tot
    52885330!
    5289 !    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( catm - ccomp_tot ) < 1.e-8 ) )  THEN
     5331!    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( conc_ijk_ugm3 - ccomp_tot ) < 1.e-8 ) )  THEN
    52905332!       !
    52915333!!--   surface concentration (almost) equal to atmospheric concentration
     
    52965338!       !
    52975339!!--    compensation point available, calculate effective resistance
    5298 !       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * catm ) / ( catm - ccomp_tot )
     5340!       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * conc_ijk_ugm3 ) / ( conc_ijk_ugm3 - ccomp_tot )
    52995341!
    53005342!    ELSE
     
    53525394 !> Boundary-layer deposition resistance following Zhang (2001)
    53535395 !------------------------------------------------------------------------
    5354  SUBROUTINE drydepo_aero_zhang_vd( vd, Rs, vs1, partsize, slipcor, nwet, tsurf, dens1, viscos1, &
    5355       luc, ftop_lu, ustar_lu )
     5396 SUBROUTINE drydepo_aero_zhang_vd( vd, rs, vs1, partsize, slipcor, nwet, tsurf, dens1, viscos1, &
     5397      luc, ftop_lu, ustar )
    53565398
    53575399    IMPLICIT NONE 
     
    53705412    REAL(wp), INTENT(IN) ::  viscos1         !< air viscosity in lowest layer
    53715413    REAL(wp), INTENT(IN) ::  ftop_lu         !< atmospheric resistnace Ra
    5372     REAL(wp), INTENT(IN) ::  ustar_lu        !< friction velocity u*   
     5414    REAL(wp), INTENT(IN) ::  ustar           !< friction velocity u*   
    53735415
    53745416    REAL(wp), INTENT(OUT) ::  vd             !< deposition velocity (m/s)
    5375     REAL(wp), INTENT(OUT) ::  Rs             !< sedimentaion resistance (s/m)
     5417    REAL(wp), INTENT(OUT) ::  rs             !< sedimentaion resistance (s/m)
    53765418!
    53775419!-- constants
     
    54175459!-- and sticking efficiency R (1 = no rebound)
    54185460    IF ( luc == ilu_ice .OR. nwet==9 .OR. luc == ilu_water_sea .OR. luc == ilu_water_inland )  THEN
    5419        stokes = vs1 * ustar_lu**2 / ( grav * kinvisc )
     5461       stokes = vs1 * ustar**2 / ( grav * kinvisc )
    54205462       Einterc = 0.0_wp
    54215463       Reffic = 1.0_wp
    54225464    ELSE IF ( luc == ilu_other .OR. luc == ilu_desert )  THEN     !<tundra of desert
    5423        stokes = vs1 * ustar_lu**2 / ( grav * kinvisc )
     5465       stokes = vs1 * ustar**2 / ( grav * kinvisc )
    54245466       Einterc = 0.0_wp
    54255467       Reffic = exp( -Stokes**0.5_wp )
    54265468    ELSE
    5427        stokes = vs1 * ustar_lu / (grav * A_lu(luc) * 1.e-3)
     5469       stokes = vs1 * ustar / (grav * A_lu(luc) * 1.e-3)
    54285470       Einterc = 0.5_wp * ( partsize / (A_lu(luc) * 1e-3 ) )**2
    54295471       Reffic = exp( -Stokes**0.5_wp )
     
    54375479!
    54385480!-- sedimentation resistance:
    5439     Rs = 1.0_wp / ( epsilon0 * ustar_lu * ( Ebrown + Eimpac + Einterc ) * Reffic )
     5481    rs = 1.0_wp / ( epsilon0 * ustar * ( Ebrown + Eimpac + Einterc ) * Reffic )
    54405482
    54415483!-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7):
     
    54475489!-- where: Rs = Rb (in Seinfeld and Pandis, 2006)
    54485490
    5449     vd = 1.0_wp / ( ftop_lu + Rs + ftop_lu * Rs * vs1) + vs1
     5491    vd = 1.0_wp / ( ftop_lu + rs + ftop_lu * rs * vs1) + vs1
    54505492
    54515493
     
    54575499 !> Original EMEP formulation by (Simpson et al, 2003) is used
    54585500 !-------------------------------------------------------------------------------------
    5459  SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb )   
     5501 SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffusivity, rb )   
    54605502
    54615503
     
    54695511    REAL(wp), INTENT(IN) ::  z0h                  !< roughness length for heat
    54705512    REAL(wp), INTENT(IN) ::  ustar                !< friction velocity
    5471     REAL(wp), INTENT(IN) ::  diffc                !< coefficient of diffusivity
    5472 
    5473     REAL(wp), INTENT(OUT) ::  Rb                  !< boundary layer resistance
     5513    REAL(wp), INTENT(IN) ::  diffusivity          !< coefficient of diffusivity
     5514
     5515    REAL(wp), INTENT(OUT) ::  rb                  !< boundary layer resistance
    54745516!
    54755517!-- const
     
    54825524!
    54835525!-- Use Simpson et al. (2003)
    5484 !-- @TODO: Check Rb over water calculation, until then leave commented lines
     5526!-- @TODO: Check rb over water calculation, until then leave commented lines
    54855527!--  IF ( is_water )  THEN
    5486 !--   org: Rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.01_wp,ustar))
    5487 !--        Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar))
     5528!--   org: rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffusivity*kappa_stab*MAX(0.01_wp,ustar))
     5529!--        rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffusivity*kappa_stab*MAX(0.1_wp,ustar))
    54885530!--  ELSE
    5489     Rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffc )**0.67_wp
     5531    rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffusivity )**0.67_wp
    54905532!--  END IF
    54915533
Note: See TracChangeset for help on using the changeset viewer.