- Timestamp:
- Apr 4, 2019 1:07:31 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/chemistry_model_mod.f90 (modified) (74 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/chemistry_model_mod.f90
r3848 r3862 27 27 ! ----------------- 28 28 ! $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 30 33 ! 31 34 ! 3784 2019-03-05 14:16:20Z banzhafs … … 2716 2719 IMPLICIT NONE 2717 2720 2718 INTEGER(iwp) :: lsp !< 2721 INTEGER(iwp) :: lsp !< running index for chem spcs. 2719 2722 2720 2723 DO lsp = 1, nspec … … 2726 2729 2727 2730 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 2728 2770 2729 2771 … … 2766 2808 INTEGER(iwp) :: k !< matching k to surface m at i,j 2767 2809 INTEGER(iwp) :: lsp !< running index for chem spcs. 2768 ! INTEGER(iwp) :: lu !< running index for landuse classes2769 2810 INTEGER(iwp) :: lu_palm !< index of PALM LSM vegetation_type at current surface element 2770 2811 INTEGER(iwp) :: lup_palm !< index of PALM LSM pavement_type at current surface element … … 2784 2825 INTEGER(iwp) :: i_pspec !< index for matching depac gas component 2785 2826 ! 2786 !-- Vegetation !<Match to DEPAC2787 INTEGER(iwp) :: ind_lu _user = 0 !< --> ERROR2788 INTEGER(iwp) :: ind_lu _b_soil = 1 !< -->ilu_desert2789 INTEGER(iwp) :: ind_lu _mixed_crops = 2 !< -->ilu_arable2790 INTEGER(iwp) :: ind_lu _s_grass = 3 !< -->ilu_grass2791 INTEGER(iwp) :: ind_lu _ev_needle_trees = 4 !< -->ilu_coniferous_forest2792 INTEGER(iwp) :: ind_lu _de_needle_trees = 5 !< -->ilu_coniferous_forest2793 INTEGER(iwp) :: ind_lu _ev_broad_trees = 6 !< -->ilu_tropical_forest2794 INTEGER(iwp) :: ind_lu _de_broad_trees = 7 !< -->ilu_deciduous_forest2795 INTEGER(iwp) :: ind_lu _t_grass = 8 !< -->ilu_grass2796 INTEGER(iwp) :: ind_lu _desert = 9 !< -->ilu_desert2797 INTEGER(iwp) :: ind_lu _tundra = 10 !< -->ilu_other2798 INTEGER(iwp) :: ind_lu _irr_crops = 11 !< -->ilu_arable2799 INTEGER(iwp) :: ind_lu _semidesert = 12 !< -->ilu_other2800 INTEGER(iwp) :: ind_lu _ice = 13 !< -->ilu_ice2801 INTEGER(iwp) :: ind_lu _marsh = 14 !< -->ilu_other2802 INTEGER(iwp) :: ind_lu _ev_shrubs = 15 !< -->ilu_mediterrean_scrub2803 INTEGER(iwp) :: ind_lu _de_shrubs = 16 !< -->ilu_mediterrean_scrub2804 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)) 2806 2847 ! 2807 2848 !-- Water 2808 INTEGER(iwp) :: ind_luw_user = 0 !< --> ERROR2809 INTEGER(iwp) :: ind_luw_lake = 1 !< -->ilu_water_inland2810 INTEGER(iwp) :: ind_luw_river = 2 !< -->ilu_water_inland2811 INTEGER(iwp) :: ind_luw_ocean = 3 !< -->ilu_water_sea2812 INTEGER(iwp) :: ind_luw_pond = 4 !< -->ilu_water_inland2813 INTEGER(iwp) :: ind_luw_fountain = 5 !< -->ilu_water_inland2849 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 2814 2855 ! 2815 2856 !-- Pavement 2816 INTEGER(iwp) :: ind_lup_user = 0 !< --> ERROR2817 INTEGER(iwp) :: ind_lup_asph_conc = 1 !< -->ilu_desert2818 INTEGER(iwp) :: ind_lup_asph = 2 !< -->ilu_desert2819 INTEGER(iwp) :: ind_lup_conc = 3 !< -->ilu_desert2820 INTEGER(iwp) :: ind_lup_sett = 4 !< -->ilu_desert2821 INTEGER(iwp) :: ind_lup_pav_stones = 5 !< -->ilu_desert2822 INTEGER(iwp) :: ind_lup_cobblest = 6 !< -->ilu_desert2823 INTEGER(iwp) :: ind_lup_metal = 7 !< -->ilu_desert2824 INTEGER(iwp) :: ind_lup_wood = 8 !< -->ilu_desert2825 INTEGER(iwp) :: ind_lup_gravel = 9 !< -->ilu_desert2826 INTEGER(iwp) :: ind_lup_f_gravel = 10 !< -->ilu_desert2827 INTEGER(iwp) :: ind_lup_pebblest = 11 !< -->ilu_desert2828 INTEGER(iwp) :: ind_lup_woodchips = 12 !< -->ilu_desert2829 INTEGER(iwp) :: ind_lup_tartan = 13 !< -->ilu_desert2830 INTEGER(iwp) :: ind_lup_art_turf = 14 !< -->ilu_desert2831 INTEGER(iwp) :: ind_lup_clay = 15 !< -->ilu_desert2857 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 2832 2873 ! 2833 2874 !-- Particle parameters according to the respective aerosol classes (PM25, PM10) … … 2836 2877 INTEGER(iwp) :: ind_p_slip = 3 !< index for slipcor in particle_pars 2837 2878 2838 INTEGER(iwp) :: part_type 2839 2840 INTEGER(iwp) :: nwet !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow2841 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 2846 2887 2847 2888 REAL(wp) :: dens !< density at layer k at i,j 2848 REAL(wp) :: r_aero_ lu!< aerodynamic resistance (s/m) at current surface element2849 REAL(wp) :: ustar_ lu!< ustar at current surface element2850 REAL(wp) :: z0h_ lu!< roughness length for heat at current surface element2851 REAL(wp) :: glrad !< rad_sw_inat current surface element2852 REAL(wp) :: ppm _to_ugm3 !< conversion factor2853 REAL(wp) :: r elh!< relative humidity at current surface element2889 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 2854 2895 REAL(wp) :: lai !< leaf area index at current surface element 2855 2896 REAL(wp) :: sai !< surface area index at current surface element assumed to be lai + 1 … … 2859 2900 REAL(wp) :: vs !< Sedimentation velocity 2860 2901 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) :: c atm !< gasconc.at i, j, k in ug/m32866 REAL(wp) :: diff c!< diffusivity2867 2868 2869 REAL(wp), DIMENSION(nspec) :: bud_ now_lu!< budget for LSM vegetation type at current surface element2870 REAL(wp), DIMENSION(nspec) :: bud_ now_lup !< budget for LSM pavement type at current surface element2871 REAL(wp), DIMENSION(nspec) :: bud_ now_luw !< budget for LSM water type at current surface element2872 REAL(wp), DIMENSION(nspec) :: bud_ now_luu !< budget for USM walls/roofs at current surface element2873 REAL(wp), DIMENSION(nspec) :: bud_ now_lug !< budget for USM green surfaces at current surface element2874 REAL(wp), DIMENSION(nspec) :: bud_ now_lud !< budget for USM windows at current surface element2875 REAL(wp), DIMENSION(nspec) :: bud _now!< overall budget at current surface element2876 REAL(wp), DIMENSION(nspec) :: c c !< concentrationi,j,k2877 REAL(wp), DIMENSION(nspec) :: ccomp_tot !< total compensation point (ug/m3)2878 !< For now kept to zero for all species!2879 2880 REAL(wp) :: t temp!< temperatur at i,j,k2902 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 2881 2922 REAL(wp) :: ts !< surface temperatur in degrees celsius 2882 2923 REAL(wp) :: qv_tmp !< surface mixing ratio at current surface element … … 3055 3096 ! 3056 3097 !-- 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) 3061 3102 lai = surf_lsm_h%lai(m) 3062 3103 sai = lai + 1 … … 3064 3105 !-- For small grid spacing neglect R_a 3065 3106 IF ( dzw(k) <= 1.0 ) THEN 3066 r_aero_ lu= 0.0_wp3107 r_aero_surf = 0.0_wp 3067 3108 ENDIF 3068 3109 ! … … 3076 3117 ! 3077 3118 !-- Initialize budgets 3078 bud_ now_lu= 0.0_wp3079 bud_ now_lup = 0.0_wp3080 bud_ now_luw = 0.0_wp3119 bud_luv = 0.0_wp 3120 bud_lup = 0.0_wp 3121 bud_luw = 0.0_wp 3081 3122 ! 3082 3123 !-- Get land use for i,j and assign to DEPAC lu 3083 3124 IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 ) THEN 3084 3125 lu_palm = surf_lsm_h%vegetation_type(m) 3085 IF ( lu_palm == ind_lu _user ) THEN3126 IF ( lu_palm == ind_luv_user ) THEN 3086 3127 message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' 3087 3128 CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 ) 3088 ELSEIF ( lu_palm == ind_lu _b_soil ) THEN3129 ELSEIF ( lu_palm == ind_luv_b_soil ) THEN 3089 3130 lu_dep = 9 3090 ELSEIF ( lu_palm == ind_lu _mixed_crops ) THEN3131 ELSEIF ( lu_palm == ind_luv_mixed_crops ) THEN 3091 3132 lu_dep = 2 3092 ELSEIF ( lu_palm == ind_lu _s_grass ) THEN3133 ELSEIF ( lu_palm == ind_luv_s_grass ) THEN 3093 3134 lu_dep = 1 3094 ELSEIF ( lu_palm == ind_lu _ev_needle_trees ) THEN3135 ELSEIF ( lu_palm == ind_luv_ev_needle_trees ) THEN 3095 3136 lu_dep = 4 3096 ELSEIF ( lu_palm == ind_lu _de_needle_trees ) THEN3137 ELSEIF ( lu_palm == ind_luv_de_needle_trees ) THEN 3097 3138 lu_dep = 4 3098 ELSEIF ( lu_palm == ind_lu _ev_broad_trees ) THEN3139 ELSEIF ( lu_palm == ind_luv_ev_broad_trees ) THEN 3099 3140 lu_dep = 12 3100 ELSEIF ( lu_palm == ind_lu _de_broad_trees ) THEN3141 ELSEIF ( lu_palm == ind_luv_de_broad_trees ) THEN 3101 3142 lu_dep = 5 3102 ELSEIF ( lu_palm == ind_lu _t_grass ) THEN3143 ELSEIF ( lu_palm == ind_luv_t_grass ) THEN 3103 3144 lu_dep = 1 3104 ELSEIF ( lu_palm == ind_lu _desert ) THEN3145 ELSEIF ( lu_palm == ind_luv_desert ) THEN 3105 3146 lu_dep = 9 3106 ELSEIF ( lu_palm == ind_lu _tundra ) THEN3147 ELSEIF ( lu_palm == ind_luv_tundra ) THEN 3107 3148 lu_dep = 8 3108 ELSEIF ( lu_palm == ind_lu _irr_crops ) THEN3149 ELSEIF ( lu_palm == ind_luv_irr_crops ) THEN 3109 3150 lu_dep = 2 3110 ELSEIF ( lu_palm == ind_lu _semidesert ) THEN3151 ELSEIF ( lu_palm == ind_luv_semidesert ) THEN 3111 3152 lu_dep = 8 3112 ELSEIF ( lu_palm == ind_lu _ice ) THEN3153 ELSEIF ( lu_palm == ind_luv_ice ) THEN 3113 3154 lu_dep = 10 3114 ELSEIF ( lu_palm == ind_lu _marsh ) THEN3155 ELSEIF ( lu_palm == ind_luv_marsh ) THEN 3115 3156 lu_dep = 8 3116 ELSEIF ( lu_palm == ind_lu _ev_shrubs ) THEN3157 ELSEIF ( lu_palm == ind_luv_ev_shrubs ) THEN 3117 3158 lu_dep = 14 3118 ELSEIF ( lu_palm == ind_lu _de_shrubs ) THEN3159 ELSEIF ( lu_palm == ind_luv_de_shrubs ) THEN 3119 3160 lu_dep = 14 3120 ELSEIF ( lu_palm == ind_lu _mixed_forest ) THEN3161 ELSEIF ( lu_palm == ind_luv_mixed_forest ) THEN 3121 3162 lu_dep = 4 3122 ELSEIF ( lu_palm == ind_lu _intrup_forest ) THEN3163 ELSEIF ( lu_palm == ind_luv_intrup_forest ) THEN 3123 3164 lu_dep = 8 3124 3165 ENDIF … … 3201 3242 !-- Concentration at i,j,k 3202 3243 DO lsp = 1, nspec 3203 c c(lsp) = chem_species(lsp)%conc(k,j,i)3244 conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i) 3204 3245 ENDDO 3205 3246 3206 3247 !-- Temperature at i,j,k 3207 t temp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp3208 3209 ts = t temp - 273.15 !< in degrees celcius3248 temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 3249 3250 ts = temp_tmp - 273.15 !< in degrees celcius 3210 3251 ! 3211 3252 !-- Viscosity of air 3212 visc = 1.496e-6 * t temp**1.5 / (ttemp+120.0)3253 visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0) 3213 3254 ! 3214 3255 !-- Air density at k … … 3217 3258 !-- Calculate relative humidity from specific humidity for DEPAC 3218 3259 qv_tmp = MAX(q(k,j,i),0.0_wp) 3219 r elh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )3260 rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) ) 3220 3261 ! 3221 3262 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget … … 3233 3274 vs = 0.0_wp 3234 3275 vd_lu = 0.0_wp 3235 Rs = 0.0_wp3236 Rb = 0.0_wp3237 Rc_tot = 0.0_wp3276 rs = 0.0_wp 3277 rb = 0.0_wp 3278 rc_tot = 0.0_wp 3238 3279 IF ( spc_names(lsp) == 'PM10' ) THEN 3239 3280 part_type = 1 … … 3245 3286 visc) 3246 3287 3247 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3288 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3248 3289 vs, & 3249 3290 particle_pars(ind_p_size, part_type), & 3250 3291 particle_pars(ind_p_slip, part_type), & 3251 nwet, t temp, dens, visc, &3292 nwet, temp_tmp, dens, visc, & 3252 3293 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) * & 3256 3297 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3257 3298 … … 3266 3307 visc) 3267 3308 3268 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3309 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3269 3310 vs, & 3270 3311 particle_pars(ind_p_size, part_type), & 3271 3312 particle_pars(ind_p_slip, part_type), & 3272 nwet, t temp, dens, visc, &3313 nwet, temp_tmp, dens, visc, & 3273 3314 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) * & 3277 3318 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3278 3319 … … 3302 3343 !-- Use density at k: 3303 3344 3304 ppm _to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m33345 ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3305 3346 ! 3306 3347 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3307 3348 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3308 c atm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m33349 conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3 3309 3350 ! 3310 3351 !-- Diffusivity for DEPAC relevant gases 3311 3352 !-- Use default value 3312 diff c= 0.11e-43353 diffusivity = 0.11e-4 3313 3354 ! 3314 3355 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3315 IF ( spc_names(lsp) == 'NO2' ) diff c= 0.136e-43316 IF ( spc_names(lsp) == 'NO' ) diff c= 0.199e-43317 IF ( spc_names(lsp) == 'O3' ) diff c= 0.144e-43318 IF ( spc_names(lsp) == 'CO' ) diff c= 0.176e-43319 IF ( spc_names(lsp) == 'SO2' ) diff c= 0.112e-43320 IF ( spc_names(lsp) == 'CH4' ) diff c= 0.191e-43321 IF ( spc_names(lsp) == 'NH3' ) diff c= 0.191e-43322 ! 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: 3324 3365 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_tot3329 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_ lu, glrad, cos_zenith, &3330 r elh, 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 ) 3332 3373 ! 3333 3374 !-- Calculate budget 3334 IF ( Rc_tot <= 0.0 ) THEN3335 3336 bud_ now_lu(lsp) = 0.0_wp3375 IF ( rc_tot <= 0.0 ) THEN 3376 3377 bud_luv(lsp) = 0.0_wp 3337 3378 3338 3379 ELSE 3339 3380 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)) * & 3343 3384 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3344 3385 ENDIF … … 3363 3404 vs = 0.0_wp 3364 3405 vd_lu = 0.0_wp 3365 Rs = 0.0_wp3366 Rb = 0.0_wp3367 Rc_tot = 0.0_wp3406 rs = 0.0_wp 3407 rb = 0.0_wp 3408 rc_tot = 0.0_wp 3368 3409 IF ( spc_names(lsp) == 'PM10' ) THEN 3369 3410 part_type = 1 … … 3375 3416 visc) 3376 3417 3377 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3418 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3378 3419 vs, & 3379 3420 particle_pars(ind_p_size, part_type), & 3380 3421 particle_pars(ind_p_slip, part_type), & 3381 nwet, t temp, dens, visc, &3422 nwet, temp_tmp, dens, visc, & 3382 3423 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) * & 3386 3427 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3387 3428 … … 3396 3437 visc) 3397 3438 3398 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3439 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3399 3440 vs, & 3400 3441 particle_pars(ind_p_size, part_type), & 3401 3442 particle_pars(ind_p_slip, part_type), & 3402 nwet, t temp, dens, visc, &3443 nwet, temp_tmp, dens, visc, & 3403 3444 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) * & 3407 3448 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3408 3449 … … 3433 3474 !-- Use density at lowest layer: 3434 3475 3435 ppm _to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m33476 ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3436 3477 ! 3437 3478 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3438 3479 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3439 c atm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m33480 conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3 3440 3481 ! 3441 3482 !-- Diffusivity for DEPAC relevant gases 3442 3483 !-- Use default value 3443 diff c= 0.11e-43484 diffusivity = 0.11e-4 3444 3485 ! 3445 3486 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3446 IF ( spc_names(lsp) == 'NO2' ) diff c= 0.136e-43447 IF ( spc_names(lsp) == 'NO' ) diff c= 0.199e-43448 IF ( spc_names(lsp) == 'O3' ) diff c= 0.144e-43449 IF ( spc_names(lsp) == 'CO' ) diff c= 0.176e-43450 IF ( spc_names(lsp) == 'SO2' ) diff c= 0.112e-43451 IF ( spc_names(lsp) == 'CH4' ) diff c= 0.191e-43452 IF ( spc_names(lsp) == 'NH3' ) diff c= 0.191e-43453 ! 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: 3455 3496 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_tot3459 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), c atm, 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 ) 3463 3504 ! 3464 3505 !-- Calculate budget 3465 IF ( Rc_tot <= 0.0 ) THEN3466 bud_ now_lup(lsp) = 0.0_wp3506 IF ( rc_tot <= 0.0 ) THEN 3507 bud_lup(lsp) = 0.0_wp 3467 3508 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)) * & 3470 3511 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3471 3512 ENDIF … … 3489 3530 vs = 0.0_wp 3490 3531 vd_lu = 0.0_wp 3491 Rs = 0.0_wp3492 Rb = 0.0_wp3493 Rc_tot = 0.0_wp3532 rs = 0.0_wp 3533 rb = 0.0_wp 3534 rc_tot = 0.0_wp 3494 3535 IF ( spc_names(lsp) == 'PM10' ) THEN 3495 3536 part_type = 1 … … 3501 3542 visc) 3502 3543 3503 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3544 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3504 3545 vs, & 3505 3546 particle_pars(ind_p_size, part_type), & 3506 3547 particle_pars(ind_p_slip, part_type), & 3507 nwet, t temp, dens, visc, &3548 nwet, temp_tmp, dens, visc, & 3508 3549 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) * & 3512 3553 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3513 3554 … … 3521 3562 visc) 3522 3563 3523 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3564 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3524 3565 vs, & 3525 3566 particle_pars(ind_p_size, part_type), & 3526 3567 particle_pars(ind_p_slip, part_type), & 3527 nwet, t temp, dens, visc, &3568 nwet, temp_tmp, dens, visc, & 3528 3569 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) * & 3532 3573 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3533 3574 … … 3558 3599 !-- Use density at lowest layer: 3559 3600 3560 ppm _to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m33601 ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3561 3602 ! 3562 3603 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3563 3604 !-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3564 c atm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m33605 conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3 3565 3606 ! 3566 3607 !-- Diffusivity for DEPAC relevant gases 3567 3608 !-- Use default value 3568 diff c= 0.11e-43609 diffusivity = 0.11e-4 3569 3610 ! 3570 3611 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3571 IF ( spc_names(lsp) == 'NO2' ) diff c= 0.136e-43572 IF ( spc_names(lsp) == 'NO' ) diff c= 0.199e-43573 IF ( spc_names(lsp) == 'O3' ) diff c= 0.144e-43574 IF ( spc_names(lsp) == 'CO' ) diff c= 0.176e-43575 IF ( spc_names(lsp) == 'SO2' ) diff c= 0.112e-43576 IF ( spc_names(lsp) == 'CH4' ) diff c= 0.191e-43577 IF ( spc_names(lsp) == 'NH3' ) diff c= 0.191e-43578 ! 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: 3580 3621 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_tot3584 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), c atm, 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 ) 3588 3629 ! 3589 3630 !-- Calculate budget 3590 IF ( Rc_tot <= 0.0 ) THEN3591 3592 bud_ now_luw(lsp) = 0.0_wp3631 IF ( rc_tot <= 0.0 ) THEN 3632 3633 bud_luw(lsp) = 0.0_wp 3593 3634 3594 3635 ELSE 3595 3636 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)) * & 3599 3640 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3600 3641 ENDIF … … 3605 3646 3606 3647 3607 bud _now= 0.0_wp3648 bud = 0.0_wp 3608 3649 ! 3609 3650 !-- Calculate overall budget for surface m and adapt concentration 3610 3651 DO lsp = 1, nspec 3611 3652 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) 3615 3656 ! 3616 3657 !-- Compute new concentration: 3617 c c(lsp) = cc(lsp) + bud_now(lsp) * inv_dh3618 3619 chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, c c(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)) 3620 3661 3621 3662 ENDDO … … 3632 3673 ! 3633 3674 !-- 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) 3638 3679 lai = surf_usm_h%lai(m) 3639 3680 sai = lai + 1 … … 3641 3682 !-- For small grid spacing neglect R_a 3642 3683 IF ( dzw(k) <= 1.0 ) THEN 3643 r_aero_ lu= 0.0_wp3684 r_aero_surf = 0.0_wp 3644 3685 ENDIF 3645 3686 ! … … 3653 3694 ! 3654 3695 !-- Initialize budgets 3655 bud_ now_luu = 0.0_wp3656 bud_ now_lug = 0.0_wp3657 bud_ now_lud = 0.0_wp3696 bud_luu = 0.0_wp 3697 bud_lug = 0.0_wp 3698 bud_lud = 0.0_wp 3658 3699 ! 3659 3700 !-- Get land use for i,j and assign to DEPAC lu … … 3662 3703 !-- For green urban surfaces (e.g. green roofs 3663 3704 !-- assume LU short grass 3664 lug_palm = ind_lu _s_grass3665 IF ( lug_palm == ind_lu _user ) THEN3705 lug_palm = ind_luv_s_grass 3706 IF ( lug_palm == ind_luv_user ) THEN 3666 3707 message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' 3667 3708 CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 ) 3668 ELSEIF ( lug_palm == ind_lu _b_soil ) THEN3709 ELSEIF ( lug_palm == ind_luv_b_soil ) THEN 3669 3710 lug_dep = 9 3670 ELSEIF ( lug_palm == ind_lu _mixed_crops ) THEN3711 ELSEIF ( lug_palm == ind_luv_mixed_crops ) THEN 3671 3712 lug_dep = 2 3672 ELSEIF ( lug_palm == ind_lu _s_grass ) THEN3713 ELSEIF ( lug_palm == ind_luv_s_grass ) THEN 3673 3714 lug_dep = 1 3674 ELSEIF ( lug_palm == ind_lu _ev_needle_trees ) THEN3715 ELSEIF ( lug_palm == ind_luv_ev_needle_trees ) THEN 3675 3716 lug_dep = 4 3676 ELSEIF ( lug_palm == ind_lu _de_needle_trees ) THEN3717 ELSEIF ( lug_palm == ind_luv_de_needle_trees ) THEN 3677 3718 lug_dep = 4 3678 ELSEIF ( lug_palm == ind_lu _ev_broad_trees ) THEN3719 ELSEIF ( lug_palm == ind_luv_ev_broad_trees ) THEN 3679 3720 lug_dep = 12 3680 ELSEIF ( lug_palm == ind_lu _de_broad_trees ) THEN3721 ELSEIF ( lug_palm == ind_luv_de_broad_trees ) THEN 3681 3722 lug_dep = 5 3682 ELSEIF ( lug_palm == ind_lu _t_grass ) THEN3723 ELSEIF ( lug_palm == ind_luv_t_grass ) THEN 3683 3724 lug_dep = 1 3684 ELSEIF ( lug_palm == ind_lu _desert ) THEN3725 ELSEIF ( lug_palm == ind_luv_desert ) THEN 3685 3726 lug_dep = 9 3686 ELSEIF ( lug_palm == ind_lu _tundra ) THEN3727 ELSEIF ( lug_palm == ind_luv_tundra ) THEN 3687 3728 lug_dep = 8 3688 ELSEIF ( lug_palm == ind_lu _irr_crops ) THEN3729 ELSEIF ( lug_palm == ind_luv_irr_crops ) THEN 3689 3730 lug_dep = 2 3690 ELSEIF ( lug_palm == ind_lu _semidesert ) THEN3731 ELSEIF ( lug_palm == ind_luv_semidesert ) THEN 3691 3732 lug_dep = 8 3692 ELSEIF ( lug_palm == ind_lu _ice ) THEN3733 ELSEIF ( lug_palm == ind_luv_ice ) THEN 3693 3734 lug_dep = 10 3694 ELSEIF ( lug_palm == ind_lu _marsh ) THEN3735 ELSEIF ( lug_palm == ind_luv_marsh ) THEN 3695 3736 lug_dep = 8 3696 ELSEIF ( lug_palm == ind_lu _ev_shrubs ) THEN3737 ELSEIF ( lug_palm == ind_luv_ev_shrubs ) THEN 3697 3738 lug_dep = 14 3698 ELSEIF ( lug_palm == ind_lu _de_shrubs ) THEN3739 ELSEIF ( lug_palm == ind_luv_de_shrubs ) THEN 3699 3740 lug_dep = 14 3700 ELSEIF ( lug_palm == ind_lu _mixed_forest ) THEN3741 ELSEIF ( lug_palm == ind_luv_mixed_forest ) THEN 3701 3742 lug_dep = 4 3702 ELSEIF ( lug_palm == ind_lu _intrup_forest ) THEN3743 ELSEIF ( lug_palm == ind_luv_intrup_forest ) THEN 3703 3744 lug_dep = 8 3704 3745 ENDIF … … 3810 3851 !-- Concentration at i,j,k 3811 3852 DO lsp = 1, nspec 3812 c c(lsp) = chem_species(lsp)%conc(k,j,i)3853 conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i) 3813 3854 ENDDO 3814 3855 ! 3815 3856 !-- Temperature at i,j,k 3816 t temp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp3817 3818 ts = t temp - 273.15 !< in degrees celcius3857 temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 3858 3859 ts = temp_tmp - 273.15 !< in degrees celcius 3819 3860 ! 3820 3861 !-- Viscosity of air 3821 visc = 1.496e-6 * t temp**1.5 / (ttemp+120.0)3862 visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0) 3822 3863 ! 3823 3864 !-- Air density at k … … 3826 3867 !-- Calculate relative humidity from specific humidity for DEPAC 3827 3868 qv_tmp = MAX(q(k,j,i),0.0_wp) 3828 r elh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )3869 rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) ) 3829 3870 ! 3830 3871 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget … … 3846 3887 vs = 0.0_wp 3847 3888 vd_lu = 0.0_wp 3848 Rs = 0.0_wp3849 Rb = 0.0_wp3850 Rc_tot = 0.0_wp3889 rs = 0.0_wp 3890 rb = 0.0_wp 3891 rc_tot = 0.0_wp 3851 3892 IF (spc_names(lsp) == 'PM10' ) THEN 3852 3893 part_type = 1 … … 3858 3899 visc) 3859 3900 3860 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3901 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3861 3902 vs, & 3862 3903 particle_pars(ind_p_size, part_type), & 3863 3904 particle_pars(ind_p_slip, part_type), & 3864 nwet, t temp, dens, visc, &3905 nwet, temp_tmp, dens, visc, & 3865 3906 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) * & 3869 3910 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3870 3911 … … 3878 3919 visc) 3879 3920 3880 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &3921 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3881 3922 vs, & 3882 3923 particle_pars(ind_p_size, part_type), & 3883 3924 particle_pars(ind_p_slip, part_type), & 3884 nwet, t temp, dens, visc, &3925 nwet, temp_tmp, dens, visc, & 3885 3926 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) * & 3889 3930 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3890 3931 … … 3914 3955 !-- Use density at k: 3915 3956 3916 ppm _to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m33957 ppm2ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3917 3958 3918 3959 ! 3919 3960 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3920 3961 !-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3921 !-- c atm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m33962 !-- conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3 3922 3963 ! 3923 3964 !-- Diffusivity for DEPAC relevant gases 3924 3965 !-- Use default value 3925 diff c= 0.11e-43966 diffusivity = 0.11e-4 3926 3967 ! 3927 3968 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3928 IF ( spc_names(lsp) == 'NO2' ) diff c= 0.136e-43929 IF ( spc_names(lsp) == 'NO' ) diff c= 0.199e-43930 IF ( spc_names(lsp) == 'O3' ) diff c= 0.144e-43931 IF ( spc_names(lsp) == 'CO' ) diff c= 0.176e-43932 IF ( spc_names(lsp) == 'SO2' ) diff c= 0.112e-43933 IF ( spc_names(lsp) == 'CH4' ) diff c= 0.191e-43934 IF ( spc_names(lsp) == 'NH3' ) diff c= 0.191e-43935 ! 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: 3937 3978 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_tot3942 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), c atm, 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 ) 3946 3987 ! 3947 3988 !-- Calculate budget 3948 IF ( Rc_tot <= 0.0 ) THEN3949 3950 bud_ now_luu(lsp) = 0.0_wp3989 IF ( rc_tot <= 0.0 ) THEN 3990 3991 bud_luu(lsp) = 0.0_wp 3951 3992 3952 3993 ELSE 3953 3994 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)) * & 3957 3998 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3958 3999 ENDIF … … 3973 4014 vs = 0.0_wp 3974 4015 vd_lu = 0.0_wp 3975 Rs = 0.0_wp3976 Rb = 0.0_wp3977 Rc_tot = 0.0_wp4016 rs = 0.0_wp 4017 rb = 0.0_wp 4018 rc_tot = 0.0_wp 3978 4019 IF ( spc_names(lsp) == 'PM10' ) THEN 3979 4020 part_type = 1 … … 3985 4026 visc) 3986 4027 3987 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &4028 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 3988 4029 vs, & 3989 4030 particle_pars(ind_p_size, part_type), & 3990 4031 particle_pars(ind_p_slip, part_type), & 3991 nwet, t temp, dens, visc, &4032 nwet, temp_tmp, dens, visc, & 3992 4033 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) * & 3996 4037 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3997 4038 … … 4005 4046 visc) 4006 4047 4007 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &4048 CALL drydepo_aero_zhang_vd( vd_lu, rs, & 4008 4049 vs, & 4009 4050 particle_pars(ind_p_size, part_type), & 4010 4051 particle_pars(ind_p_slip, part_type), & 4011 nwet, t temp, dens, visc, &4052 nwet, temp_tmp, dens, visc, & 4012 4053 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) * & 4016 4057 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4017 4058 … … 4041 4082 !-- Use density at k: 4042 4083 4043 ppm _to_ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m34084 ppm2ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m3 4044 4085 ! 4045 4086 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4046 4087 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4047 c atm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m34088 conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3 4048 4089 ! 4049 4090 !-- Diffusivity for DEPAC relevant gases 4050 4091 !-- Use default value 4051 diff c= 0.11e-44092 diffusivity = 0.11e-4 4052 4093 ! 4053 4094 !-- overwrite with known coefficients of diffusivity from Massman (1998) 4054 IF ( spc_names(lsp) == 'NO2' ) diff c= 0.136e-44055 IF ( spc_names(lsp) == 'NO' ) diff c= 0.199e-44056 IF ( spc_names(lsp) == 'O3' ) diff c= 0.144e-44057 IF ( spc_names(lsp) == 'CO' ) diff c= 0.176e-44058 IF ( spc_names(lsp) == 'SO2' ) diff c= 0.112e-44059 IF ( spc_names(lsp) == 'CH4' ) diff c= 0.191e-44060 IF ( spc_names(lsp) == 'NH3' ) diff c= 0.191e-44061 ! 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: 4063 4104 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_tot4068 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), c atm, 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 ) 4072 4113 ! 4073 4114 !-- Calculate budget 4074 IF ( Rc_tot <= 0.0 ) THEN4075 4076 bud_ now_lug(lsp) = 0.0_wp4115 IF ( rc_tot <= 0.0 ) THEN 4116 4117 bud_lug(lsp) = 0.0_wp 4077 4118 4078 4119 ELSE 4079 4120 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)) * & 4083 4124 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4084 4125 ENDIF … … 4103 4144 vs = 0.0_wp 4104 4145 vd_lu = 0.0_wp 4105 Rs = 0.0_wp4106 Rb = 0.0_wp4107 Rc_tot = 0.0_wp4146 rs = 0.0_wp 4147 rb = 0.0_wp 4148 rc_tot = 0.0_wp 4108 4149 IF ( spc_names(lsp) == 'PM10' ) THEN 4109 4150 part_type = 1 … … 4118 4159 particle_pars(ind_p_size, part_type), & 4119 4160 particle_pars(ind_p_slip, part_type), & 4120 nwet, t temp, 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) * & 4124 4165 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4125 4166 … … 4136 4177 particle_pars(ind_p_size, part_type), & 4137 4178 particle_pars(ind_p_slip, part_type), & 4138 nwet, t temp, dens, visc, &4179 nwet, temp_tmp, dens, visc, & 4139 4180 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) * & 4143 4184 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4144 4185 … … 4168 4209 !-- Use density at k: 4169 4210 4170 ppm _to_ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m34211 ppm2ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m3 4171 4212 ! 4172 4213 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4173 4214 !-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4174 c atm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m34215 conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 * specmolm(i_pspec) ! in ug/m3 4175 4216 ! 4176 4217 !-- Diffusivity for DEPAC relevant gases 4177 4218 !-- Use default value 4178 diff c= 0.11e-44219 diffusivity = 0.11e-4 4179 4220 ! 4180 4221 !-- overwrite with known coefficients of diffusivity from Massman (1998) 4181 IF ( spc_names(lsp) == 'NO2' ) diff c= 0.136e-44182 IF ( spc_names(lsp) == 'NO' ) diff c= 0.199e-44183 IF ( spc_names(lsp) == 'O3' ) diff c= 0.144e-44184 IF ( spc_names(lsp) == 'CO' ) diff c= 0.176e-44185 IF ( spc_names(lsp) == 'SO2' ) diff c= 0.112e-44186 IF ( spc_names(lsp) == 'CH4' ) diff c= 0.191e-44187 IF ( spc_names(lsp) == 'NH3' ) diff c= 0.191e-44188 ! 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: 4190 4231 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_tot4194 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), c atm, 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 ) 4198 4239 ! 4199 4240 !-- Calculate budget 4200 IF ( Rc_tot <= 0.0 ) THEN4201 4202 bud_ now_lud(lsp) = 0.0_wp4241 IF ( rc_tot <= 0.0 ) THEN 4242 4243 bud_lud(lsp) = 0.0_wp 4203 4244 4204 4245 ELSE 4205 4246 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)) * & 4209 4250 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4210 4251 ENDIF … … 4215 4256 4216 4257 4217 bud _now= 0.0_wp4258 bud = 0.0_wp 4218 4259 ! 4219 4260 !-- Calculate overall budget for surface m and adapt concentration … … 4221 4262 4222 4263 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) 4226 4267 ! 4227 4268 !-- Compute new concentration 4228 c c(lsp) = cc(lsp) + bud_now(lsp) * inv_dh4229 4230 chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, c c(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) ) 4231 4272 4232 4273 ENDDO … … 4250 4291 !> van Zanten et al., 2010. 4251 4292 !------------------------------------------------------------------------------! 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, c atm, 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, & 4254 4295 ra, rb ) 4255 4296 ! … … 4291 4332 REAL(wp), INTENT(IN) :: t !< temperature (C) 4292 4333 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) 4294 4335 REAL(wp), INTENT(IN) :: sinphi !< sin of solar elevation angle 4295 4336 REAL(wp), INTENT(IN) :: rh !< relative humidity (%) 4296 4337 REAL(wp), INTENT(IN) :: lai !< one-sidedleaf area index (-) 4297 4338 REAL(wp), INTENT(IN) :: sai !< surface area index (-) (lai + branches and stems) 4298 REAL(wp), INTENT(IN) :: diff c!< diffusivity4339 REAL(wp), INTENT(IN) :: diffusivity !< diffusivity 4299 4340 REAL(wp), INTENT(IN) :: p !< pressure (Pa) 4300 REAL(wp), INTENT(IN) :: c atm !< actual atmospheric concentration (ug/m3)4341 REAL(wp), INTENT(IN) :: conc_ijk_ugm3 !< actual atmospheric concentration (ug/m3), in DEPAC=Catm 4301 4342 REAL(wp), INTENT(IN) :: ra !< aerodynamic resistance (s/m) 4302 4343 REAL(wp), INTENT(IN) :: rb !< boundary layer resistance (s/m) … … 4340 4381 ! 4341 4382 !-- Next statement is just to avoid compiler warning about unused variable 4342 IF ( day_of_year == 0 .OR. ( c atm+ lat + ra + rb ) > 0.0_wp ) CONTINUE4383 IF ( day_of_year == 0 .OR. ( conc_ijk_ugm3 + lat + ra + rb ) > 0.0_wp ) CONTINUE 4343 4384 ! 4344 4385 !-- Define component number … … 4406 4447 ! 4407 4448 !-- 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 ) 4409 4450 ! 4410 4451 !-- Effective soil conductance: … … 4419 4460 !-- lai_present, sai_present, & 4420 4461 !-- ccomp_tot, & 4421 !-- c atm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &4462 !-- conc_ijk_ugm3=conc_ijk_ugm3,c_ave_prev_nh3=c_ave_prev_nh3, & 4422 4463 !-- c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, & 4423 4464 !-- tsea=tsea,cw=cw,cstom=cstom,csoil=csoil ) … … 4426 4467 !-- IF ( present(rc_eff) ) then 4427 4468 !-- check on required arguments: 4428 !-- IF ( (.not. present(c atm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then4429 !-- stop 'output argument rc_eff requires input arguments c atm, 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' 4430 4471 !-- END IF 4431 4472 ! 4432 4473 !-- Compute rc_eff : 4433 ! CALL rc_comp_point_rc_eff(ccomp_tot,c atm,ra,rb,rc_tot,rc_eff)4474 ! CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff) 4434 4475 ! ENDIF 4435 4476 ENDIF … … 4714 4755 !> Subroutine to compute stomatal conductance 4715 4756 !------------------------------------------------------------------------------! 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 ) 4717 4758 4718 4759 IMPLICIT NONE … … 4727 4768 4728 4769 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) 4730 4771 REAL(wp), INTENT(IN) :: sinphi !< sin of solar elevation angle 4731 4772 REAL(wp), INTENT(IN) :: t !< temperature (C) 4732 4773 REAL(wp), INTENT(IN) :: rh !< relative humidity (%) 4733 REAL(wp), INTENT(IN) :: diff c!< diffusion coefficient of the gas involved4774 REAL(wp), INTENT(IN) :: diffusivity !< diffusion coefficient of the gas involved 4734 4775 4735 4776 REAL(wp), OPTIONAL,INTENT(IN) :: p !< pressure (Pa) … … 4757 4798 IF ( lai_present ) THEN 4758 4799 4759 IF ( glrad > 0.0_wp ) THEN4800 IF ( solar_rad > 0.0_wp ) THEN 4760 4801 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 * diff c/ dO3 !< Gstom of Emberson is derived for ozone4802 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 4763 4804 ELSE 4764 4805 gstom = 0.0_wp … … 4783 4824 !> Subroutine to compute stomatal conductance according to Emberson 4784 4825 !------------------------------------------------------------------------------! 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 ) 4786 4827 ! 4787 4828 !> History … … 4823 4864 LOGICAL, INTENT(IN) :: lai_present 4824 4865 4825 REAL(wp), INTENT(IN) :: glrad !< global radiation(W/m2)4866 REAL(wp), INTENT(IN) :: solar_rad !< solar radiation, dirict+diffuse (W/m2) 4826 4867 REAL(wp), INTENT(IN) :: t !< temperature (C) 4827 4868 REAL(wp), INTENT(IN) :: vpd !< vapour pressure deficit (kPa) … … 4872 4913 ! 4873 4914 !-- 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 ) 4875 4916 ! 4876 4917 !-- par for shaded leaves (canopy averaged): 4877 4918 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 ) THEN4919 IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp ) THEN 4879 4920 parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi ) !< Zhang et al., 2001 4880 4921 END IF … … 4883 4924 !-- alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5 4884 4925 parsun = pardir * 0.5/sinphi + parshade !< Norman, 1982 4885 IF ( glrad > 200.0_wp .AND. lai > 2.5_wp ) THEN4926 IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp ) THEN 4886 4927 parsun = pardir**0.8 * 0.5 / sinphi + parshade !< Zhang et al., 2001 4887 4928 END IF … … 4945 4986 !> Meteorological Service of Canada 4946 4987 !> 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) 4948 4989 !> 4949 4990 !> @todo Check/connect/replace with radiation_model_mod variables 4950 4991 !------------------------------------------------------------------- 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 ) 4952 4993 4953 4994 IMPLICIT NONE 4954 4995 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) 4956 4997 REAL(wp), INTENT(IN) :: sinphi !< sine of the solar elevation 4957 4998 REAL(wp), INTENT(IN) :: pres !< actual pressure (to correct for height) (Pa) … … 4995 5036 rn = MAX( 0.01_wp, rdm + rdn ) 4996 5037 ! 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 ) ) 4999 5041 ! 5000 5042 !-- Calculate total visible radiation … … 5015 5057 !> rc_get_vpd: get vapour pressure deficit (kPa) 5016 5058 !------------------------------------------------------------------- 5017 SUBROUTINE rc_get_vpd( temp, r elh, vpd )5059 SUBROUTINE rc_get_vpd( temp, rh, vpd ) 5018 5060 5019 5061 IMPLICIT NONE … … 5021 5063 !-- Input/output variables: 5022 5064 REAL(wp), INTENT(IN) :: temp !< temperature (C) 5023 REAL(wp), INTENT(IN) :: r elh !< relative humidity (%)5065 REAL(wp), INTENT(IN) :: rh !< relative humidity (%) 5024 5066 5025 5067 REAL(wp), INTENT(OUT) :: vpd !< vapour pressure deficit (kPa) … … 5038 5080 !-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973) 5039 5081 esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5 5040 vpd = esat * ( 1 - r elh / 100 )5082 vpd = esat * ( 1 - rh / 100 ) 5041 5083 5042 5084 END SUBROUTINE rc_get_vpd … … 5267 5309 !> 5268 5310 ! ------------------------------------------------------------------------------------------- 5269 ! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, c atm, ra, rb, rc_tot, rc_eff )5311 ! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, conc_ijk_ugm3, ra, rb, rc_tot, rc_eff ) 5270 5312 ! 5271 5313 ! IMPLICIT NONE 5272 5314 ! 5273 5315 !!-- 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) :: c atm !< 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) 5281 5323 ! 5282 5324 ! ! … … 5287 5329 ! rc_eff = rc_tot 5288 5330 ! 5289 ! ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( c atm- ccomp_tot ) < 1.e-8 ) ) THEN5331 ! ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( conc_ijk_ugm3 - ccomp_tot ) < 1.e-8 ) ) THEN 5290 5332 ! ! 5291 5333 !!-- surface concentration (almost) equal to atmospheric concentration … … 5296 5338 ! ! 5297 5339 !!-- compensation point available, calculate effective resistance 5298 ! rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * c atm ) / ( catm- ccomp_tot )5340 ! rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * conc_ijk_ugm3 ) / ( conc_ijk_ugm3 - ccomp_tot ) 5299 5341 ! 5300 5342 ! ELSE … … 5352 5394 !> Boundary-layer deposition resistance following Zhang (2001) 5353 5395 !------------------------------------------------------------------------ 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 ) 5356 5398 5357 5399 IMPLICIT NONE … … 5370 5412 REAL(wp), INTENT(IN) :: viscos1 !< air viscosity in lowest layer 5371 5413 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* 5373 5415 5374 5416 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) 5376 5418 ! 5377 5419 !-- constants … … 5417 5459 !-- and sticking efficiency R (1 = no rebound) 5418 5460 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 ) 5420 5462 Einterc = 0.0_wp 5421 5463 Reffic = 1.0_wp 5422 5464 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 ) 5424 5466 Einterc = 0.0_wp 5425 5467 Reffic = exp( -Stokes**0.5_wp ) 5426 5468 ELSE 5427 stokes = vs1 * ustar _lu/ (grav * A_lu(luc) * 1.e-3)5469 stokes = vs1 * ustar / (grav * A_lu(luc) * 1.e-3) 5428 5470 Einterc = 0.5_wp * ( partsize / (A_lu(luc) * 1e-3 ) )**2 5429 5471 Reffic = exp( -Stokes**0.5_wp ) … … 5437 5479 ! 5438 5480 !-- 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 ) 5440 5482 5441 5483 !-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7): … … 5447 5489 !-- where: Rs = Rb (in Seinfeld and Pandis, 2006) 5448 5490 5449 vd = 1.0_wp / ( ftop_lu + Rs + ftop_lu * Rs * vs1) + vs15491 vd = 1.0_wp / ( ftop_lu + rs + ftop_lu * rs * vs1) + vs1 5450 5492 5451 5493 … … 5457 5499 !> Original EMEP formulation by (Simpson et al, 2003) is used 5458 5500 !------------------------------------------------------------------------------------- 5459 SUBROUTINE get_rb_cell( is_water, z0h, ustar, diff c, Rb )5501 SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffusivity, rb ) 5460 5502 5461 5503 … … 5469 5511 REAL(wp), INTENT(IN) :: z0h !< roughness length for heat 5470 5512 REAL(wp), INTENT(IN) :: ustar !< friction velocity 5471 REAL(wp), INTENT(IN) :: diff c!< coefficient of diffusivity5472 5473 REAL(wp), INTENT(OUT) :: Rb !< boundary layer resistance5513 REAL(wp), INTENT(IN) :: diffusivity !< coefficient of diffusivity 5514 5515 REAL(wp), INTENT(OUT) :: rb !< boundary layer resistance 5474 5516 ! 5475 5517 !-- const … … 5482 5524 ! 5483 5525 !-- Use Simpson et al. (2003) 5484 !-- @TODO: Check Rb over water calculation, until then leave commented lines5526 !-- @TODO: Check rb over water calculation, until then leave commented lines 5485 5527 !-- 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)) 5488 5530 !-- ELSE 5489 Rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffc)**0.67_wp5531 rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffusivity )**0.67_wp 5490 5532 !-- END IF 5491 5533
Note: See TracChangeset
for help on using the changeset viewer.
