| 2570 | |
| 2571 | |
| 2572 | |
| 2573 | !-------------------------------------------------------------------------------! |
| 2574 | ! |
| 2575 | ! Description: |
| 2576 | ! ------------ |
| 2577 | !> Subroutine to calculate the deposition of gases and PMs. For now deposition |
| 2578 | !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT |
| 2579 | !> considered. The deposition of particles is derived following Zhang et al., |
| 2580 | !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010). |
| 2581 | !> |
| 2582 | !> @TODO: Consider deposition on vertical surfaces |
| 2583 | !> @TODO: Consider overlaying horizontal surfaces |
| 2584 | !> @TODO: Check error messages |
| 2585 | !-------------------------------------------------------------------------------! |
| 2586 | |
| 2587 | SUBROUTINE chem_depo (i,j) |
| 2588 | |
| 2589 | |
| 2590 | USE control_parameters, & |
| 2591 | ONLY: dt_3d, intermediate_timestep_count, latitude |
| 2592 | |
| 2593 | USE arrays_3d, & |
| 2594 | ONLY: dzw, rho_air_zw |
| 2595 | |
| 2596 | USE date_and_time_mod, & |
| 2597 | ONLY: day_of_year |
| 2598 | |
| 2599 | USE surface_mod, & |
| 2600 | ONLY: surf_lsm_h, surf_usm_h, surf_type, ind_veg_wall, & |
| 2601 | ind_pav_green, ind_wat_win |
| 2602 | |
| 2603 | USE radiation_model_mod, & |
| 2604 | ONLY: zenith |
| 2605 | |
| 2606 | |
| 2607 | |
| 2608 | IMPLICIT NONE |
| 2609 | |
| 2610 | INTEGER,INTENT(IN) :: i,j |
| 2611 | |
| 2612 | INTEGER(iwp) :: k ! matching k to surface m at i,j |
| 2613 | INTEGER(iwp) :: lsp ! running index for chem spcs. |
| 2614 | INTEGER(iwp) :: lu ! running index for landuse classes |
| 2615 | INTEGER(iwp) :: lu_palm ! index of PALM LSM vegetation_type at current surface element |
| 2616 | INTEGER(iwp) :: lup_palm ! index of PALM LSM pavement_type at current surface element |
| 2617 | INTEGER(iwp) :: luw_palm ! index of PALM LSM water_type at current surface element |
| 2618 | INTEGER(iwp) :: luu_palm ! index of PALM USM walls/roofs at current surface element |
| 2619 | INTEGER(iwp) :: lug_palm ! index of PALM USM green walls/roofs at current surface element |
| 2620 | INTEGER(iwp) :: lud_palm ! index of PALM USM windows at current surface element |
| 2621 | INTEGER(iwp) :: lu_dep ! matching DEPAC LU to lu_palm |
| 2622 | INTEGER(iwp) :: lup_dep ! matching DEPAC LU to lup_palm |
| 2623 | INTEGER(iwp) :: luw_dep ! matching DEPAC LU to luw_palm |
| 2624 | INTEGER(iwp) :: luu_dep ! matching DEPAC LU to luu_palm |
| 2625 | INTEGER(iwp) :: lug_dep ! matching DEPAC LU to lug_palm |
| 2626 | INTEGER(iwp) :: lud_dep ! matching DEPAC LU to lud_palm |
| 2627 | INTEGER(iwp) :: m ! index for horizontal surfaces |
| 2628 | |
| 2629 | INTEGER(iwp) :: pspec ! running index |
| 2630 | INTEGER(iwp) :: i_pspec ! index for matching depac gas component |
| 2631 | |
| 2632 | |
| 2633 | ! Vegetation !Match to DEPAC |
| 2634 | INTEGER(iwp) :: ind_lu_user = 0 ! --> ERROR |
| 2635 | INTEGER(iwp) :: ind_lu_b_soil = 1 ! --> ilu_desert |
| 2636 | INTEGER(iwp) :: ind_lu_mixed_crops = 2 ! --> ilu_arable |
| 2637 | INTEGER(iwp) :: ind_lu_s_grass = 3 ! --> ilu_grass |
| 2638 | INTEGER(iwp) :: ind_lu_ev_needle_trees = 4 ! --> ilu_coniferous_forest |
| 2639 | INTEGER(iwp) :: ind_lu_de_needle_trees = 5 ! --> ilu_coniferous_forest |
| 2640 | INTEGER(iwp) :: ind_lu_ev_broad_trees = 6 ! --> ilu_tropical_forest |
| 2641 | INTEGER(iwp) :: ind_lu_de_broad_trees = 7 ! --> ilu_deciduous_forest |
| 2642 | INTEGER(iwp) :: ind_lu_t_grass = 8 ! --> ilu_grass |
| 2643 | INTEGER(iwp) :: ind_lu_desert = 9 ! --> ilu_desert |
| 2644 | INTEGER(iwp) :: ind_lu_tundra = 10 ! --> ilu_other |
| 2645 | INTEGER(iwp) :: ind_lu_irr_crops = 11 ! --> ilu_arable |
| 2646 | INTEGER(iwp) :: ind_lu_semidesert = 12 ! --> ilu_other |
| 2647 | INTEGER(iwp) :: ind_lu_ice = 13 ! --> ilu_ice |
| 2648 | INTEGER(iwp) :: ind_lu_marsh = 14 ! --> ilu_other |
| 2649 | INTEGER(iwp) :: ind_lu_ev_shrubs = 15 ! --> ilu_mediterrean_scrub |
| 2650 | INTEGER(iwp) :: ind_lu_de_shrubs = 16 ! --> ilu_mediterrean_scrub |
| 2651 | INTEGER(iwp) :: ind_lu_mixed_forest = 17 ! --> ilu_coniferous_forest (ave(decid+conif)) |
| 2652 | INTEGER(iwp) :: ind_lu_intrup_forest = 18 ! --> ilu_other (ave(other+decid)) |
| 2653 | |
| 2654 | ! Water |
| 2655 | INTEGER(iwp) :: ind_luw_user = 0 ! --> ERROR |
| 2656 | INTEGER(iwp) :: ind_luw_lake = 1 ! --> ilu_water_inland |
| 2657 | INTEGER(iwp) :: ind_luw_river = 2 ! --> ilu_water_inland |
| 2658 | INTEGER(iwp) :: ind_luw_ocean = 3 ! --> ilu_water_sea |
| 2659 | INTEGER(iwp) :: ind_luw_pond = 4 ! --> ilu_water_inland |
| 2660 | INTEGER(iwp) :: ind_luw_fountain = 5 ! --> ilu_water_inland |
| 2661 | |
| 2662 | ! Pavement |
| 2663 | INTEGER(iwp) :: ind_lup_user = 0 ! --> ERROR |
| 2664 | INTEGER(iwp) :: ind_lup_asph_conc = 1 ! --> ilu_desert |
| 2665 | INTEGER(iwp) :: ind_lup_asph = 2 ! --> ilu_desert |
| 2666 | INTEGER(iwp) :: ind_lup_conc = 3 ! --> ilu_desert |
| 2667 | INTEGER(iwp) :: ind_lup_sett = 4 ! --> ilu_desert |
| 2668 | INTEGER(iwp) :: ind_lup_pav_stones = 5 ! --> ilu_desert |
| 2669 | INTEGER(iwp) :: ind_lup_cobblest = 6 ! --> ilu_desert |
| 2670 | INTEGER(iwp) :: ind_lup_metal = 7 ! --> ilu_desert |
| 2671 | INTEGER(iwp) :: ind_lup_wood = 8 ! --> ilu_desert |
| 2672 | INTEGER(iwp) :: ind_lup_gravel = 9 ! --> ilu_desert |
| 2673 | INTEGER(iwp) :: ind_lup_f_gravel = 10 ! --> ilu_desert |
| 2674 | INTEGER(iwp) :: ind_lup_pebblest = 11 ! --> ilu_desert |
| 2675 | INTEGER(iwp) :: ind_lup_woodchips = 12 ! --> ilu_desert |
| 2676 | INTEGER(iwp) :: ind_lup_tartan = 13 ! --> ilu_desert |
| 2677 | INTEGER(iwp) :: ind_lup_art_turf = 14 ! --> ilu_desert |
| 2678 | INTEGER(iwp) :: ind_lup_clay = 15 ! --> ilu_desert |
| 2679 | |
| 2680 | |
| 2681 | |
| 2682 | !-- Particle parameters according to the respective aerosol classes (PM25, PM10) |
| 2683 | INTEGER(iwp) :: ind_p_size = 1 !< index for partsize in particle_pars |
| 2684 | INTEGER(iwp) :: ind_p_dens = 2 !< index for rhopart in particle_pars |
| 2685 | INTEGER(iwp) :: ind_p_slip = 3 !< index for slipcor in particle_pars |
| 2686 | |
| 2687 | INTEGER(iwp) :: part_type |
| 2688 | |
| 2689 | INTEGER(iwp) :: nwet ! wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow |
| 2690 | |
| 2691 | REAL(kind=wp) :: dt_chem |
| 2692 | REAL(kind=wp) :: dh |
| 2693 | REAL(kind=wp) :: inv_dh |
| 2694 | REAL(kind=wp) :: dt_dh |
| 2695 | |
| 2696 | REAL(kind=wp) :: dens ! Density at lowest layer at i,j |
| 2697 | REAL(kind=wp) :: r_aero_lu ! aerodynamic resistance (s/m) at current surface element |
| 2698 | REAL(kind=wp) :: ustar_lu ! ustar at current surface element |
| 2699 | REAL(kind=wp) :: z0h_lu ! roughness length for heat at current surface element |
| 2700 | REAL(kind=wp) :: glrad ! rad_sw_in at current surface element |
| 2701 | REAL(kind=wp) :: ppm_to_ugm3 ! conversion factor |
| 2702 | REAL(kind=wp) :: relh ! relative humidity at current surface element |
| 2703 | REAL(kind=wp) :: lai ! leaf area index at current surface element |
| 2704 | REAL(kind=wp) :: sai ! surface area index at current surface element assumed to be lai + 1 |
| 2705 | |
| 2706 | REAL(kind=wp) :: slinnfac |
| 2707 | REAL(kind=wp) :: visc ! Viscosity |
| 2708 | REAL(kind=wp) :: vs ! Sedimentation velocity |
| 2709 | REAL(kind=wp) :: vd_lu ! deposition velocity (m/s) |
| 2710 | REAL(kind=wp) :: Rs ! Sedimentaion resistance (s/m) |
| 2711 | REAL(kind=wp) :: Rb ! quasi-laminar boundary layer resistance (s/m) |
| 2712 | REAL(kind=wp) :: Rc_tot ! total canopy resistance Rc (s/m) |
| 2713 | |
| 2714 | REAL(kind=wp) :: catm ! gasconc. at i,j,k=1 in ug/m3 |
| 2715 | REAL(kind=wp) :: diffc ! Diffusivity |
| 2716 | |
| 2717 | |
| 2718 | REAL(kind=wp),DIMENSION(nspec) :: bud_now_lu ! budget for LSM vegetation type at current surface element |
| 2719 | REAL(kind=wp),DIMENSION(nspec) :: bud_now_lup ! budget for LSM pavement type at current surface element |
| 2720 | REAL(kind=wp),DIMENSION(nspec) :: bud_now_luw ! budget for LSM water type at current surface element |
| 2721 | REAL(kind=wp),DIMENSION(nspec) :: bud_now_luu ! budget for USM walls/roofs at current surface element |
| 2722 | REAL(kind=wp),DIMENSION(nspec) :: bud_now_lug ! budget for USM green surfaces at current surface element |
| 2723 | REAL(kind=wp),DIMENSION(nspec) :: bud_now_lud ! budget for USM windows at current surface element |
| 2724 | REAL(kind=wp),DIMENSION(nspec) :: bud_now ! overall budget at current surface element |
| 2725 | REAL(kind=wp),DIMENSION(nspec) :: cc ! concentration i,j,k |
| 2726 | REAL(kind=wp),DIMENSION(nspec) :: ccomp_tot ! total compensation point (ug/m3) (= 0 for species that don't have a compensation) |
| 2727 | ! For now kept to zero for all species! |
| 2728 | |
| 2729 | REAL(kind=wp) :: ttemp ! temperatur at i,j,k |
| 2730 | REAL(kind=wp) :: ts ! surface temperatur in degrees celsius |
| 2731 | REAL(kind=wp) :: qv_tmp ! surface mixing ratio at current surface element |
| 2732 | |
| 2733 | |
| 2734 | !-- Particle parameters (PM10 (1), PM25 (2)) |
| 2735 | !-- partsize (diameter in m) , rhopart (density in kg/m3), slipcor (slip correction factor DIMENSIONless, Seinfeld and Pandis 2006, Table 9.3) |
| 2736 | REAL(wp), DIMENSION(1:3,1:2), PARAMETER :: particle_pars = RESHAPE( (/ & |
| 2737 | 8.0e-6_wp, 1.14e3_wp, 1.016_wp, & ! 1 |
| 2738 | 0.7e-6_wp, 1.14e3_wp, 1.082_wp & ! 2 |
| 2739 | /), (/ 3, 2 /) ) |
| 2740 | |
| 2741 | |
| 2742 | LOGICAL :: match_lsm ! flag indicating natural-type surface |
| 2743 | LOGICAL :: match_usm ! flag indicating urban-type surface |
| 2744 | |
| 2745 | |
| 2746 | |
| 2747 | !-- List of names of possible tracers |
| 2748 | CHARACTER(len=*), PARAMETER :: pspecnames(69) = (/ & |
| 2749 | 'NO2 ', & ! NO2 |
| 2750 | 'NO ', & ! NO |
| 2751 | 'O3 ', & ! O3 |
| 2752 | 'CO ', & ! CO |
| 2753 | 'form ', & ! FORM |
| 2754 | 'ald ', & ! ALD |
| 2755 | 'pan ', & ! PAN |
| 2756 | 'mgly ', & ! MGLY |
| 2757 | 'par ', & ! PAR |
| 2758 | 'ole ', & ! OLE |
| 2759 | 'eth ', & ! ETH |
| 2760 | 'tol ', & ! TOL |
| 2761 | 'cres ', & ! CRES |
| 2762 | 'xyl ', & ! XYL |
| 2763 | 'SO4a_f ', & ! SO4a_f |
| 2764 | 'SO2 ', & ! SO2 |
| 2765 | 'HNO2 ', & ! HNO2 |
| 2766 | 'CH4 ', & ! CH4 |
| 2767 | 'NH3 ', & ! NH3 |
| 2768 | 'NO3 ', & ! NO3 |
| 2769 | 'OH ', & ! OH |
| 2770 | 'HO2 ', & ! HO2 |
| 2771 | 'N2O5 ', & ! N2O5 |
| 2772 | 'SO4a_c ', & ! SO4a_c |
| 2773 | 'NH4a_f ', & ! NH4a_f |
| 2774 | 'NO3a_f ', & ! NO3a_f |
| 2775 | 'NO3a_c ', & ! NO3a_c |
| 2776 | 'C2O3 ', & ! C2O3 |
| 2777 | 'XO2 ', & ! XO2 |
| 2778 | 'XO2N ', & ! XO2N |
| 2779 | 'cro ', & ! CRO |
| 2780 | 'HNO3 ', & ! HNO3 |
| 2781 | 'H2O2 ', & ! H2O2 |
| 2782 | 'iso ', & ! ISO |
| 2783 | 'ispd ', & ! ISPD |
| 2784 | 'to2 ', & ! TO2 |
| 2785 | 'open ', & ! OPEN |
| 2786 | 'terp ', & ! TERP |
| 2787 | 'ec_f ', & ! EC_f |
| 2788 | 'ec_c ', & ! EC_c |
| 2789 | 'pom_f ', & ! POM_f |
| 2790 | 'pom_c ', & ! POM_c |
| 2791 | 'ppm_f ', & ! PPM_f |
| 2792 | 'ppm_c ', & ! PPM_c |
| 2793 | 'na_ff ', & ! Na_ff |
| 2794 | 'na_f ', & ! Na_f |
| 2795 | 'na_c ', & ! Na_c |
| 2796 | 'na_cc ', & ! Na_cc |
| 2797 | 'na_ccc ', & ! Na_ccc |
| 2798 | 'dust_ff ', & ! dust_ff |
| 2799 | 'dust_f ', & ! dust_f |
| 2800 | 'dust_c ', & ! dust_c |
| 2801 | 'dust_cc ', & ! dust_cc |
| 2802 | 'dust_ccc ', & ! dust_ccc |
| 2803 | 'tpm10 ', & ! tpm10 |
| 2804 | 'tpm25 ', & ! tpm25 |
| 2805 | 'tss ', & ! tss |
| 2806 | 'tdust ', & ! tdust |
| 2807 | 'tc ', & ! tc |
| 2808 | 'tcg ', & ! tcg |
| 2809 | 'tsoa ', & ! tsoa |
| 2810 | 'tnmvoc ', & ! tnmvoc |
| 2811 | 'SOxa ', & ! SOxa |
| 2812 | 'NOya ', & ! NOya |
| 2813 | 'NHxa ', & ! NHxa |
| 2814 | 'NO2_obs ', & ! NO2_obs |
| 2815 | 'tpm10_biascorr', & ! tpm10_biascorr |
| 2816 | 'tpm25_biascorr', & ! tpm25_biascorr |
| 2817 | 'O3_biascorr ' /) ! o3_biascorr |
| 2818 | |
| 2819 | |
| 2820 | |
| 2821 | ! tracer mole mass: |
| 2822 | REAL, PARAMETER :: specmolm(69) = (/ & |
| 2823 | xm_O * 2 + xm_N, & ! NO2 |
| 2824 | xm_O + xm_N, & ! NO |
| 2825 | xm_O * 3, & ! O3 |
| 2826 | xm_C + xm_O, & ! CO |
| 2827 | xm_H * 2 + xm_C + xm_O, & ! FORM |
| 2828 | xm_H * 3 + xm_C * 2 + xm_O, & ! ALD |
| 2829 | xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, & ! PAN |
| 2830 | xm_H * 4 + xm_C * 3 + xm_O * 2, & ! MGLY |
| 2831 | xm_H * 3 + xm_C, & ! PAR |
| 2832 | xm_H * 3 + xm_C * 2, & ! OLE |
| 2833 | xm_H * 4 + xm_C * 2, & ! ETH |
| 2834 | xm_H * 8 + xm_C * 7, & ! TOL |
| 2835 | xm_H * 8 + xm_C * 7 + xm_O, & ! CRES |
| 2836 | xm_H * 10 + xm_C * 8, & ! XYL |
| 2837 | xm_S + xm_O * 4, & ! SO4a_f |
| 2838 | xm_S + xm_O * 2, & ! SO2 |
| 2839 | xm_H + xm_O * 2 + xm_N, & ! HNO2 |
| 2840 | xm_H * 4 + xm_C, & ! CH4 |
| 2841 | xm_H * 3 + xm_N, & ! NH3 |
| 2842 | xm_O * 3 + xm_N, & ! NO3 |
| 2843 | xm_H + xm_O, & ! OH |
| 2844 | xm_H + xm_O * 2, & ! HO2 |
| 2845 | xm_O * 5 + xm_N * 2, & ! N2O5 |
| 2846 | xm_S + xm_O * 4, & ! SO4a_c |
| 2847 | xm_H * 4 + xm_N, & ! NH4a_f |
| 2848 | xm_O * 3 + xm_N, & ! NO3a_f |
| 2849 | xm_O * 3 + xm_N, & ! NO3a_c |
| 2850 | xm_C * 2 + xm_O * 3, & ! C2O3 |
| 2851 | xm_dummy, & ! XO2 |
| 2852 | xm_dummy, & ! XO2N |
| 2853 | xm_dummy, & ! CRO |
| 2854 | xm_H + xm_O * 3 + xm_N, & ! HNO3 |
| 2855 | xm_H * 2 + xm_O * 2, & ! H2O2 |
| 2856 | xm_H * 8 + xm_C * 5, & ! ISO |
| 2857 | xm_dummy, & ! ISPD |
| 2858 | xm_dummy, & ! TO2 |
| 2859 | xm_dummy, & ! OPEN |
| 2860 | xm_H * 16 + xm_C * 10, & ! TERP |
| 2861 | xm_dummy, & ! EC_f |
| 2862 | xm_dummy, & ! EC_c |
| 2863 | xm_dummy, & ! POM_f |
| 2864 | xm_dummy, & ! POM_c |
| 2865 | xm_dummy, & ! PPM_f |
| 2866 | xm_dummy, & ! PPM_c |
| 2867 | xm_Na, & ! Na_ff |
| 2868 | xm_Na, & ! Na_f |
| 2869 | xm_Na, & ! Na_c |
| 2870 | xm_Na, & ! Na_cc |
| 2871 | xm_Na, & ! Na_ccc |
| 2872 | xm_dummy, & ! dust_ff |
| 2873 | xm_dummy, & ! dust_f |
| 2874 | xm_dummy, & ! dust_c |
| 2875 | xm_dummy, & ! dust_cc |
| 2876 | xm_dummy, & ! dust_ccc |
| 2877 | xm_dummy, & ! tpm10 |
| 2878 | xm_dummy, & ! tpm25 |
| 2879 | xm_dummy, & ! tss |
| 2880 | xm_dummy, & ! tdust |
| 2881 | xm_dummy, & ! tc |
| 2882 | xm_dummy, & ! tcg |
| 2883 | xm_dummy, & ! tsoa |
| 2884 | xm_dummy, & ! tnmvoc |
| 2885 | xm_dummy, & ! SOxa |
| 2886 | xm_dummy, & ! NOya |
| 2887 | xm_dummy, & ! NHxa |
| 2888 | xm_O * 2 + xm_N, & ! NO2_obs |
| 2889 | xm_dummy, & ! tpm10_biascorr |
| 2890 | xm_dummy, & ! tpm25_biascorr |
| 2891 | xm_O * 3 /) ! o3_biascorr |
| 2892 | |
| 2893 | |
| 2894 | |
| 2895 | ! Initialize surface element m |
| 2896 | m = 0 |
| 2897 | k = 0 |
| 2898 | ! LSM or USM surface present at i,j: |
| 2899 | ! Default surfaces are NOT considered for deposition |
| 2900 | match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) |
| 2901 | match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) |
| 2902 | |
| 2903 | |
| 2904 | !!!!!!! |
| 2905 | ! LSM ! |
| 2906 | !!!!!!! |
| 2907 | |
| 2908 | IF ( match_lsm ) THEN |
| 2909 | |
| 2910 | ! Get surface element information at i,j: |
| 2911 | m = surf_lsm_h%start_index(j,i) |
| 2912 | |
| 2913 | ! Get corresponding k |
| 2914 | k = surf_lsm_h%k(m) |
| 2915 | |
| 2916 | ! Get needed variables for surface element m |
| 2917 | ustar_lu = surf_lsm_h%us(m) |
| 2918 | z0h_lu = surf_lsm_h%z0h(m) |
| 2919 | r_aero_lu = surf_lsm_h%r_a(m) |
| 2920 | glrad = surf_lsm_h%rad_sw_in(m) |
| 2921 | lai = surf_lsm_h%lai(m) |
| 2922 | sai = lai + 1 |
| 2923 | |
| 2924 | ! For small grid spacing neglect R_a |
| 2925 | IF (dzw(k) <= 1.0) THEN |
| 2926 | r_aero_lu = 0.0_wp |
| 2927 | ENDIF |
| 2928 | |
| 2929 | !Initialize lu's |
| 2930 | lu_palm = 0 |
| 2931 | lu_dep = 0 |
| 2932 | lup_palm = 0 |
| 2933 | lup_dep = 0 |
| 2934 | luw_palm = 0 |
| 2935 | luw_dep = 0 |
| 2936 | |
| 2937 | !Initialize budgets |
| 2938 | bud_now_lu = 0.0_wp |
| 2939 | bud_now_lup = 0.0_wp |
| 2940 | bud_now_luw = 0.0_wp |
| 2941 | |
| 2942 | |
| 2943 | ! Get land use for i,j and assign to DEPAC lu |
| 2944 | IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN |
| 2945 | lu_palm = surf_lsm_h%vegetation_type(m) |
| 2946 | IF (lu_palm == ind_lu_user) THEN |
| 2947 | message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' |
| 2948 | CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 ) |
| 2949 | ELSEIF (lu_palm == ind_lu_b_soil) THEN |
| 2950 | lu_dep = 9 |
| 2951 | ELSEIF (lu_palm == ind_lu_mixed_crops) THEN |
| 2952 | lu_dep = 2 |
| 2953 | ELSEIF (lu_palm == ind_lu_s_grass) THEN |
| 2954 | lu_dep = 1 |
| 2955 | ELSEIF (lu_palm == ind_lu_ev_needle_trees) THEN |
| 2956 | lu_dep = 4 |
| 2957 | ELSEIF (lu_palm == ind_lu_de_needle_trees) THEN |
| 2958 | lu_dep = 4 |
| 2959 | ELSEIF (lu_palm == ind_lu_ev_broad_trees) THEN |
| 2960 | lu_dep = 12 |
| 2961 | ELSEIF (lu_palm == ind_lu_de_broad_trees) THEN |
| 2962 | lu_dep = 5 |
| 2963 | ELSEIF (lu_palm == ind_lu_t_grass) THEN |
| 2964 | lu_dep = 1 |
| 2965 | ELSEIF (lu_palm == ind_lu_desert) THEN |
| 2966 | lu_dep = 9 |
| 2967 | ELSEIF (lu_palm == ind_lu_tundra ) THEN |
| 2968 | lu_dep = 8 |
| 2969 | ELSEIF (lu_palm == ind_lu_irr_crops) THEN |
| 2970 | lu_dep = 2 |
| 2971 | ELSEIF (lu_palm == ind_lu_semidesert) THEN |
| 2972 | lu_dep = 8 |
| 2973 | ELSEIF (lu_palm == ind_lu_ice) THEN |
| 2974 | lu_dep = 10 |
| 2975 | ELSEIF (lu_palm == ind_lu_marsh) THEN |
| 2976 | lu_dep = 8 |
| 2977 | ELSEIF (lu_palm == ind_lu_ev_shrubs) THEN |
| 2978 | lu_dep = 14 |
| 2979 | ELSEIF (lu_palm == ind_lu_de_shrubs ) THEN |
| 2980 | lu_dep = 14 |
| 2981 | ELSEIF (lu_palm == ind_lu_mixed_forest) THEN |
| 2982 | lu_dep = 4 |
| 2983 | ELSEIF (lu_palm == ind_lu_intrup_forest) THEN |
| 2984 | lu_dep = 8 |
| 2985 | ENDIF |
| 2986 | ENDIF |
| 2987 | |
| 2988 | IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN |
| 2989 | lup_palm = surf_lsm_h%pavement_type(m) |
| 2990 | IF (lup_palm == ind_lup_user) THEN |
| 2991 | message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' |
| 2992 | CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 ) |
| 2993 | ELSEIF (lup_palm == ind_lup_asph_conc) THEN |
| 2994 | lup_dep = 9 |
| 2995 | ELSEIF (lup_palm == ind_lup_asph) THEN |
| 2996 | lup_dep = 9 |
| 2997 | ELSEIF (lup_palm == ind_lup_conc) THEN |
| 2998 | lup_dep = 9 |
| 2999 | ELSEIF (lup_palm == ind_lup_sett) THEN |
| 3000 | lup_dep = 9 |
| 3001 | ELSEIF (lup_palm == ind_lup_pav_stones) THEN |
| 3002 | lup_dep = 9 |
| 3003 | ELSEIF (lup_palm == ind_lup_cobblest) THEN |
| 3004 | lup_dep = 9 |
| 3005 | ELSEIF (lup_palm == ind_lup_metal) THEN |
| 3006 | lup_dep = 9 |
| 3007 | ELSEIF (lup_palm == ind_lup_wood) THEN |
| 3008 | lup_dep = 9 |
| 3009 | ELSEIF (lup_palm == ind_lup_gravel) THEN |
| 3010 | lup_dep = 9 |
| 3011 | ELSEIF (lup_palm == ind_lup_f_gravel) THEN |
| 3012 | lup_dep = 9 |
| 3013 | ELSEIF (lup_palm == ind_lup_pebblest) THEN |
| 3014 | lup_dep = 9 |
| 3015 | ELSEIF (lup_palm == ind_lup_woodchips) THEN |
| 3016 | lup_dep = 9 |
| 3017 | ELSEIF (lup_palm == ind_lup_tartan) THEN |
| 3018 | lup_dep = 9 |
| 3019 | ELSEIF (lup_palm == ind_lup_art_turf) THEN |
| 3020 | lup_dep = 9 |
| 3021 | ELSEIF (lup_palm == ind_lup_clay) THEN |
| 3022 | lup_dep = 9 |
| 3023 | ENDIF |
| 3024 | ENDIF |
| 3025 | |
| 3026 | IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN |
| 3027 | luw_palm = surf_lsm_h%water_type(m) |
| 3028 | IF (luw_palm == ind_luw_user) THEN |
| 3029 | message_string = 'No water type defined. Please define water type to enable deposition calculation' |
| 3030 | CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 ) |
| 3031 | ELSEIF (luw_palm == ind_luw_lake) THEN |
| 3032 | luw_dep = 13 |
| 3033 | ELSEIF (luw_palm == ind_luw_river) THEN |
| 3034 | luw_dep = 13 |
| 3035 | ELSEIF (luw_palm == ind_luw_ocean) THEN |
| 3036 | luw_dep = 6 |
| 3037 | ELSEIF (luw_palm == ind_luw_pond) THEN |
| 3038 | luw_dep = 13 |
| 3039 | ELSEIF (luw_palm == ind_luw_fountain) THEN |
| 3040 | luw_dep = 13 |
| 3041 | ENDIF |
| 3042 | ENDIF |
| 3043 | |
| 3044 | |
| 3045 | |
| 3046 | ! Set wetness indicator to dry or wet for lsm vegetation or pavement |
| 3047 | IF (surf_lsm_h%c_liq(m) > 0) THEN |
| 3048 | nwet = 1 |
| 3049 | ELSE |
| 3050 | nwet = 0 |
| 3051 | ENDIF |
| 3052 | |
| 3053 | ! Compute length of time step |
| 3054 | IF ( call_chem_at_all_substeps ) THEN |
| 3055 | dt_chem = dt_3d * weight_pres(intermediate_timestep_count) |
| 3056 | ELSE |
| 3057 | dt_chem = dt_3d |
| 3058 | ENDIF |
| 3059 | |
| 3060 | |
| 3061 | dh = dzw(k) |
| 3062 | inv_dh = 1.0_wp / dh |
| 3063 | dt_dh = dt_chem/dh |
| 3064 | |
| 3065 | ! Concentration at i,j,k |
| 3066 | DO lsp = 1, nspec |
| 3067 | cc(lsp) = chem_species(lsp)%conc(k,j,i) |
| 3068 | ENDDO |
| 3069 | |
| 3070 | |
| 3071 | ! Temperature column at i,j |
| 3072 | ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp |
| 3073 | |
| 3074 | ! Surface temperature in degrees Celcius: |
| 3075 | ts = ttemp - 273.15 ! in degrees celcius |
| 3076 | |
| 3077 | ! Viscosity of air in lowest layer |
| 3078 | visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) |
| 3079 | |
| 3080 | ! Air density in lowest layer |
| 3081 | dens = rho_air_zw(k) |
| 3082 | |
| 3083 | ! Calculate relative humidity from specific humidity for DEPAC |
| 3084 | qv_tmp = max(q(k,j,i),0.0_wp) |
| 3085 | relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) ) |
| 3086 | |
| 3087 | |
| 3088 | |
| 3089 | ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget |
| 3090 | ! for each surface fraction. Then derive overall budget taking into account the surface fractions. |
| 3091 | |
| 3092 | ! Vegetation |
| 3093 | IF (surf_lsm_h%frac(ind_veg_wall,m) > 0) THEN |
| 3094 | |
| 3095 | |
| 3096 | slinnfac = 1.0_wp |
| 3097 | |
| 3098 | ! Get vd |
| 3099 | DO lsp = 1, nvar |
| 3100 | !Initialize |
| 3101 | vs = 0.0_wp |
| 3102 | vd_lu = 0.0_wp |
| 3103 | Rs = 0.0_wp |
| 3104 | Rb = 0.0_wp |
| 3105 | Rc_tot = 0.0_wp |
| 3106 | IF (spc_names(lsp) == 'PM10') THEN |
| 3107 | part_type = 1 |
| 3108 | ! sedimentation velocity in lowest layer |
| 3109 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3110 | particle_pars(ind_p_size, part_type), & |
| 3111 | particle_pars(ind_p_slip, part_type), & |
| 3112 | visc) |
| 3113 | |
| 3114 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3115 | vs, & |
| 3116 | particle_pars(ind_p_size, part_type), & |
| 3117 | particle_pars(ind_p_slip, part_type), & |
| 3118 | nwet, ttemp, dens, visc, & |
| 3119 | lu_dep, & |
| 3120 | r_aero_lu, ustar_lu) |
| 3121 | |
| 3122 | bud_now_lu(lsp) = - cc(lsp) * & |
| 3123 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3124 | |
| 3125 | |
| 3126 | ELSEIF (spc_names(lsp) == 'PM25') THEN |
| 3127 | part_type = 2 |
| 3128 | ! sedimentation velocity in lowest layer: |
| 3129 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3130 | particle_pars(ind_p_size, part_type), & |
| 3131 | particle_pars(ind_p_slip, part_type), & |
| 3132 | visc) |
| 3133 | |
| 3134 | |
| 3135 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3136 | vs, & |
| 3137 | particle_pars(ind_p_size, part_type), & |
| 3138 | particle_pars(ind_p_slip, part_type), & |
| 3139 | nwet, ttemp, dens, visc, & |
| 3140 | lu_dep , & |
| 3141 | r_aero_lu, ustar_lu) |
| 3142 | |
| 3143 | bud_now_lu(lsp) = - cc(lsp) * & |
| 3144 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3145 | |
| 3146 | |
| 3147 | ELSE |
| 3148 | |
| 3149 | ! GASES |
| 3150 | |
| 3151 | ! Read spc_name of current species for gas parameter |
| 3152 | |
| 3153 | IF (ANY(pspecnames(:) == spc_names(lsp))) THEN |
| 3154 | i_pspec = 0 |
| 3155 | DO pspec = 1, 69 |
| 3156 | IF (pspecnames(pspec) == spc_names(lsp)) THEN |
| 3157 | i_pspec = pspec |
| 3158 | END IF |
| 3159 | ENDDO |
| 3160 | |
| 3161 | ELSE |
| 3162 | ! Species for now not deposited |
| 3163 | CYCLE |
| 3164 | ENDIF |
| 3165 | |
| 3166 | ! Factor used for conversion from ppb to ug/m3 : |
| 3167 | ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & |
| 3168 | ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) |
| 3169 | ! c 1e-9 xm_tracer 1e9 / xm_air dens |
| 3170 | ! thus: |
| 3171 | ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 |
| 3172 | ! Use density at lowest layer: |
| 3173 | |
| 3174 | ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 |
| 3175 | |
| 3176 | ! atmospheric concentration in DEPAC is requested in ug/m3: |
| 3177 | ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole |
| 3178 | catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 |
| 3179 | |
| 3180 | ! Diffusivity for DEPAC relevant gases |
| 3181 | ! Use default value |
| 3182 | diffc = 0.11e-4 |
| 3183 | ! overwrite with known coefficients of diffusivity from Massman (1998) |
| 3184 | IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 |
| 3185 | IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 |
| 3186 | IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 |
| 3187 | IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 |
| 3188 | IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 |
| 3189 | IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 |
| 3190 | IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 |
| 3191 | |
| 3192 | |
| 3193 | ! Get quasi-laminar boundary layer resistance Rb: |
| 3194 | CALL get_rb_cell( (lu_dep == ilu_water_sea) .or. (lu_dep == ilu_water_inland), & |
| 3195 | z0h_lu, ustar_lu, diffc, & |
| 3196 | Rb) |
| 3197 | |
| 3198 | ! Get Rc_tot |
| 3199 | CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & |
| 3200 | relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & |
| 3201 | r_aero_lu , Rb) |
| 3202 | |
| 3203 | |
| 3204 | ! Calculate budget |
| 3205 | IF (Rc_tot .le. 0.0) THEN |
| 3206 | |
| 3207 | bud_now_lu(lsp) = 0.0_wp |
| 3208 | |
| 3209 | ELSE |
| 3210 | |
| 3211 | ! Compute exchange velocity for current lu: |
| 3212 | vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) |
| 3213 | |
| 3214 | bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & |
| 3215 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3216 | ENDIF |
| 3217 | |
| 3218 | ENDIF |
| 3219 | ENDDO |
| 3220 | ENDIF |
| 3221 | |
| 3222 | |
| 3223 | ! Pavement |
| 3224 | IF (surf_lsm_h%frac(ind_pav_green,m) > 0) THEN |
| 3225 | |
| 3226 | |
| 3227 | ! No vegetation on pavements: |
| 3228 | lai = 0.0_wp |
| 3229 | sai = 0.0_wp |
| 3230 | |
| 3231 | slinnfac = 1.0_wp |
| 3232 | |
| 3233 | ! Get vd |
| 3234 | DO lsp = 1, nvar |
| 3235 | !Initialize |
| 3236 | vs = 0.0_wp |
| 3237 | vd_lu = 0.0_wp |
| 3238 | Rs = 0.0_wp |
| 3239 | Rb = 0.0_wp |
| 3240 | Rc_tot = 0.0_wp |
| 3241 | IF (spc_names(lsp) == 'PM10') THEN |
| 3242 | part_type = 1 |
| 3243 | ! sedimentation velocity in lowest layer: |
| 3244 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3245 | particle_pars(ind_p_size, part_type), & |
| 3246 | particle_pars(ind_p_slip, part_type), & |
| 3247 | visc) |
| 3248 | |
| 3249 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3250 | vs, & |
| 3251 | particle_pars(ind_p_size, part_type), & |
| 3252 | particle_pars(ind_p_slip, part_type), & |
| 3253 | nwet, ttemp, dens, visc, & |
| 3254 | lup_dep, & |
| 3255 | r_aero_lu, ustar_lu) |
| 3256 | |
| 3257 | bud_now_lup(lsp) = - cc(lsp) * & |
| 3258 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3259 | |
| 3260 | |
| 3261 | ELSEIF (spc_names(lsp) == 'PM25') THEN |
| 3262 | part_type = 2 |
| 3263 | ! sedimentation velocity in lowest layer: |
| 3264 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3265 | particle_pars(ind_p_size, part_type), & |
| 3266 | particle_pars(ind_p_slip, part_type), & |
| 3267 | visc) |
| 3268 | |
| 3269 | |
| 3270 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3271 | vs, & |
| 3272 | particle_pars(ind_p_size, part_type), & |
| 3273 | particle_pars(ind_p_slip, part_type), & |
| 3274 | nwet, ttemp, dens, visc, & |
| 3275 | lup_dep, & |
| 3276 | r_aero_lu, ustar_lu) |
| 3277 | |
| 3278 | bud_now_lup(lsp) = - cc(lsp) * & |
| 3279 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3280 | |
| 3281 | |
| 3282 | ELSE |
| 3283 | |
| 3284 | ! GASES |
| 3285 | |
| 3286 | ! Read spc_name of current species for gas parameter |
| 3287 | |
| 3288 | IF (ANY(pspecnames(:) == spc_names(lsp))) THEN |
| 3289 | i_pspec = 0 |
| 3290 | DO pspec = 1, 69 |
| 3291 | IF (pspecnames(pspec) == spc_names(lsp)) THEN |
| 3292 | i_pspec = pspec |
| 3293 | END IF |
| 3294 | ENDDO |
| 3295 | |
| 3296 | ELSE |
| 3297 | ! Species for now not deposited |
| 3298 | CYCLE |
| 3299 | ENDIF |
| 3300 | |
| 3301 | ! Factor used for conversion from ppb to ug/m3 : |
| 3302 | ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & |
| 3303 | ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) |
| 3304 | ! c 1e-9 xm_tracer 1e9 / xm_air dens |
| 3305 | ! thus: |
| 3306 | ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 |
| 3307 | ! Use density at lowest layer: |
| 3308 | |
| 3309 | ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 |
| 3310 | |
| 3311 | ! atmospheric concentration in DEPAC is requested in ug/m3: |
| 3312 | ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole |
| 3313 | catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 |
| 3314 | |
| 3315 | ! Diffusivity for DEPAC relevant gases |
| 3316 | ! Use default value |
| 3317 | diffc = 0.11e-4 |
| 3318 | ! overwrite with known coefficients of diffusivity from Massman (1998) |
| 3319 | IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 |
| 3320 | IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 |
| 3321 | IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 |
| 3322 | IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 |
| 3323 | IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 |
| 3324 | IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 |
| 3325 | IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 |
| 3326 | |
| 3327 | |
| 3328 | ! Get quasi-laminar boundary layer resistance Rb: |
| 3329 | CALL get_rb_cell( (lup_dep == ilu_water_sea) .or. (lup_dep == ilu_water_inland), & |
| 3330 | z0h_lu, ustar_lu, diffc, & |
| 3331 | Rb) |
| 3332 | |
| 3333 | |
| 3334 | ! Get Rc_tot |
| 3335 | CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & |
| 3336 | relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & |
| 3337 | r_aero_lu , Rb) |
| 3338 | |
| 3339 | |
| 3340 | ! Calculate budget |
| 3341 | IF (Rc_tot .le. 0.0) THEN |
| 3342 | |
| 3343 | bud_now_lup(lsp) = 0.0_wp |
| 3344 | |
| 3345 | ELSE |
| 3346 | |
| 3347 | ! Compute exchange velocity for current lu: |
| 3348 | vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) |
| 3349 | |
| 3350 | bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & |
| 3351 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3352 | ENDIF |
| 3353 | |
| 3354 | |
| 3355 | ENDIF |
| 3356 | ENDDO |
| 3357 | ENDIF |
| 3358 | |
| 3359 | |
| 3360 | ! Water |
| 3361 | IF (surf_lsm_h%frac(ind_wat_win,m) > 0) THEN |
| 3362 | |
| 3363 | |
| 3364 | ! No vegetation on water: |
| 3365 | lai = 0.0_wp |
| 3366 | sai = 0.0_wp |
| 3367 | |
| 3368 | slinnfac = 1.0_wp |
| 3369 | |
| 3370 | ! Get vd |
| 3371 | DO lsp = 1, nvar |
| 3372 | !Initialize |
| 3373 | vs = 0.0_wp |
| 3374 | vd_lu = 0.0_wp |
| 3375 | Rs = 0.0_wp |
| 3376 | Rb = 0.0_wp |
| 3377 | Rc_tot = 0.0_wp |
| 3378 | IF (spc_names(lsp) == 'PM10') THEN |
| 3379 | part_type = 1 |
| 3380 | ! sedimentation velocity in lowest layer: |
| 3381 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3382 | particle_pars(ind_p_size, part_type), & |
| 3383 | particle_pars(ind_p_slip, part_type), & |
| 3384 | visc) |
| 3385 | |
| 3386 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3387 | vs, & |
| 3388 | particle_pars(ind_p_size, part_type), & |
| 3389 | particle_pars(ind_p_slip, part_type), & |
| 3390 | nwet, ttemp, dens, visc, & |
| 3391 | luw_dep, & |
| 3392 | r_aero_lu, ustar_lu) |
| 3393 | |
| 3394 | bud_now_luw(lsp) = - cc(lsp) * & |
| 3395 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3396 | |
| 3397 | |
| 3398 | ELSEIF (spc_names(lsp) == 'PM25') THEN |
| 3399 | part_type = 2 |
| 3400 | ! sedimentation velocity in lowest layer: |
| 3401 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3402 | particle_pars(ind_p_size, part_type), & |
| 3403 | particle_pars(ind_p_slip, part_type), & |
| 3404 | visc) |
| 3405 | |
| 3406 | |
| 3407 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3408 | vs, & |
| 3409 | particle_pars(ind_p_size, part_type), & |
| 3410 | particle_pars(ind_p_slip, part_type), & |
| 3411 | nwet, ttemp, dens, visc, & |
| 3412 | luw_dep, & |
| 3413 | r_aero_lu, ustar_lu) |
| 3414 | |
| 3415 | bud_now_luw(lsp) = - cc(lsp) * & |
| 3416 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3417 | |
| 3418 | |
| 3419 | ELSE |
| 3420 | |
| 3421 | ! GASES |
| 3422 | |
| 3423 | ! Read spc_name of current species for gas PARAMETER |
| 3424 | |
| 3425 | IF (ANY(pspecnames(:) == spc_names(lsp))) THEN |
| 3426 | i_pspec = 0 |
| 3427 | DO pspec = 1, 69 |
| 3428 | IF (pspecnames(pspec) == spc_names(lsp)) THEN |
| 3429 | i_pspec = pspec |
| 3430 | END IF |
| 3431 | ENDDO |
| 3432 | |
| 3433 | ELSE |
| 3434 | ! Species for now not deposited |
| 3435 | CYCLE |
| 3436 | ENDIF |
| 3437 | |
| 3438 | ! Factor used for conversion from ppb to ug/m3 : |
| 3439 | ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & |
| 3440 | ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) |
| 3441 | ! c 1e-9 xm_tracer 1e9 / xm_air dens |
| 3442 | ! thus: |
| 3443 | ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 |
| 3444 | ! Use density at lowest layer: |
| 3445 | |
| 3446 | ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 |
| 3447 | |
| 3448 | ! atmospheric concentration in DEPAC is requested in ug/m3: |
| 3449 | ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole |
| 3450 | catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 |
| 3451 | |
| 3452 | ! Diffusivity for DEPAC relevant gases |
| 3453 | ! Use default value |
| 3454 | diffc = 0.11e-4 |
| 3455 | ! overwrite with known coefficients of diffusivity from Massman (1998) |
| 3456 | IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 |
| 3457 | IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 |
| 3458 | IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 |
| 3459 | IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 |
| 3460 | IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 |
| 3461 | IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 |
| 3462 | IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 |
| 3463 | |
| 3464 | |
| 3465 | ! Get quasi-laminar boundary layer resistance Rb: |
| 3466 | CALL get_rb_cell( (luw_dep == ilu_water_sea) .or. (luw_dep == ilu_water_inland), & |
| 3467 | z0h_lu, ustar_lu, diffc, & |
| 3468 | Rb) |
| 3469 | |
| 3470 | ! Get Rc_tot |
| 3471 | CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & |
| 3472 | relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & |
| 3473 | r_aero_lu , Rb) |
| 3474 | |
| 3475 | |
| 3476 | ! Calculate budget |
| 3477 | IF (Rc_tot .le. 0.0) THEN |
| 3478 | |
| 3479 | bud_now_luw(lsp) = 0.0_wp |
| 3480 | |
| 3481 | ELSE |
| 3482 | |
| 3483 | ! Compute exchange velocity for current lu: |
| 3484 | vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) |
| 3485 | |
| 3486 | bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & |
| 3487 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3488 | ENDIF |
| 3489 | |
| 3490 | ENDIF |
| 3491 | ENDDO |
| 3492 | ENDIF |
| 3493 | |
| 3494 | |
| 3495 | bud_now = 0.0_wp |
| 3496 | ! Calculate budget for surface m and adapt concentration |
| 3497 | DO lsp = 1, nspec |
| 3498 | |
| 3499 | |
| 3500 | bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m)*bud_now_lu(lsp) + & |
| 3501 | surf_lsm_h%frac(ind_pav_green,m)*bud_now_lup(lsp) + & |
| 3502 | surf_lsm_h%frac(ind_wat_win,m)*bud_now_luw(lsp) |
| 3503 | |
| 3504 | ! Compute new concentration, add contribution of all exchange processes: |
| 3505 | cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh |
| 3506 | |
| 3507 | chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp)) |
| 3508 | |
| 3509 | ENDDO |
| 3510 | |
| 3511 | ENDIF |
| 3512 | |
| 3513 | |
| 3514 | |
| 3515 | !!!!!!! |
| 3516 | ! USM ! |
| 3517 | !!!!!!! |
| 3518 | |
| 3519 | IF ( match_usm ) THEN |
| 3520 | |
| 3521 | ! Get surface element information at i,j: |
| 3522 | m = surf_usm_h%start_index(j,i) |
| 3523 | |
| 3524 | k = surf_usm_h%k(m) |
| 3525 | |
| 3526 | ! Get needed variables for surface element m |
| 3527 | ustar_lu = surf_usm_h%us(m) |
| 3528 | z0h_lu = surf_usm_h%z0h(m) |
| 3529 | r_aero_lu = surf_usm_h%r_a(m) |
| 3530 | glrad = surf_usm_h%rad_sw_in(m) |
| 3531 | lai = surf_usm_h%lai(m) |
| 3532 | sai = lai + 1 |
| 3533 | |
| 3534 | ! For small grid spacing neglect R_a |
| 3535 | IF (dzw(k) <= 1.0) THEN |
| 3536 | r_aero_lu = 0.0_wp |
| 3537 | ENDIF |
| 3538 | |
| 3539 | !Initialize lu's |
| 3540 | luu_palm = 0 |
| 3541 | luu_dep = 0 |
| 3542 | lug_palm = 0 |
| 3543 | lug_dep = 0 |
| 3544 | lud_palm = 0 |
| 3545 | lud_dep = 0 |
| 3546 | |
| 3547 | !Initialize budgets |
| 3548 | bud_now_luu = 0.0_wp |
| 3549 | bud_now_lug = 0.0_wp |
| 3550 | bud_now_lud = 0.0_wp |
| 3551 | |
| 3552 | |
| 3553 | ! Get land use for i,j and assign to DEPAC lu |
| 3554 | IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN |
| 3555 | ! For green urban surfaces (e.g. green roofs |
| 3556 | ! assume LU short grass |
| 3557 | lug_palm = ind_lu_s_grass |
| 3558 | IF (lug_palm == ind_lu_user) THEN |
| 3559 | message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation' |
| 3560 | CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 ) |
| 3561 | ELSEIF (lug_palm == ind_lu_b_soil) THEN |
| 3562 | lug_dep = 9 |
| 3563 | ELSEIF (lug_palm == ind_lu_mixed_crops) THEN |
| 3564 | lug_dep = 2 |
| 3565 | ELSEIF (lug_palm == ind_lu_s_grass) THEN |
| 3566 | lug_dep = 1 |
| 3567 | ELSEIF (lug_palm == ind_lu_ev_needle_trees) THEN |
| 3568 | lug_dep = 4 |
| 3569 | ELSEIF (lug_palm == ind_lu_de_needle_trees) THEN |
| 3570 | lug_dep = 4 |
| 3571 | ELSEIF (lug_palm == ind_lu_ev_broad_trees) THEN |
| 3572 | lug_dep = 12 |
| 3573 | ELSEIF (lug_palm == ind_lu_de_broad_trees) THEN |
| 3574 | lug_dep = 5 |
| 3575 | ELSEIF (lug_palm == ind_lu_t_grass) THEN |
| 3576 | lug_dep = 1 |
| 3577 | ELSEIF (lug_palm == ind_lu_desert) THEN |
| 3578 | lug_dep = 9 |
| 3579 | ELSEIF (lug_palm == ind_lu_tundra ) THEN |
| 3580 | lug_dep = 8 |
| 3581 | ELSEIF (lug_palm == ind_lu_irr_crops) THEN |
| 3582 | lug_dep = 2 |
| 3583 | ELSEIF (lug_palm == ind_lu_semidesert) THEN |
| 3584 | lug_dep = 8 |
| 3585 | ELSEIF (lug_palm == ind_lu_ice) THEN |
| 3586 | lug_dep = 10 |
| 3587 | ELSEIF (lug_palm == ind_lu_marsh) THEN |
| 3588 | lug_dep = 8 |
| 3589 | ELSEIF (lug_palm == ind_lu_ev_shrubs) THEN |
| 3590 | lug_dep = 14 |
| 3591 | ELSEIF (lug_palm == ind_lu_de_shrubs ) THEN |
| 3592 | lug_dep = 14 |
| 3593 | ELSEIF (lug_palm == ind_lu_mixed_forest) THEN |
| 3594 | lug_dep = 4 |
| 3595 | ELSEIF (lug_palm == ind_lu_intrup_forest) THEN |
| 3596 | lug_dep = 8 |
| 3597 | ENDIF |
| 3598 | ENDIF |
| 3599 | |
| 3600 | IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN |
| 3601 | ! For walls in USM assume concrete walls/roofs, |
| 3602 | ! assumed LU class desert as also assumed for |
| 3603 | ! pavements in LSM |
| 3604 | luu_palm = ind_lup_conc |
| 3605 | IF (luu_palm == ind_lup_user) THEN |
| 3606 | message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' |
| 3607 | CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 ) |
| 3608 | ELSEIF (luu_palm == ind_lup_asph_conc) THEN |
| 3609 | luu_dep = 9 |
| 3610 | ELSEIF (luu_palm == ind_lup_asph) THEN |
| 3611 | luu_dep = 9 |
| 3612 | ELSEIF (luu_palm == ind_lup_conc) THEN |
| 3613 | luu_dep = 9 |
| 3614 | ELSEIF (luu_palm == ind_lup_sett) THEN |
| 3615 | luu_dep = 9 |
| 3616 | ELSEIF (luu_palm == ind_lup_pav_stones) THEN |
| 3617 | luu_dep = 9 |
| 3618 | ELSEIF (luu_palm == ind_lup_cobblest) THEN |
| 3619 | luu_dep = 9 |
| 3620 | ELSEIF (luu_palm == ind_lup_metal) THEN |
| 3621 | luu_dep = 9 |
| 3622 | ELSEIF (luu_palm == ind_lup_wood) THEN |
| 3623 | luu_dep = 9 |
| 3624 | ELSEIF (luu_palm == ind_lup_gravel) THEN |
| 3625 | luu_dep = 9 |
| 3626 | ELSEIF (luu_palm == ind_lup_f_gravel) THEN |
| 3627 | luu_dep = 9 |
| 3628 | ELSEIF (luu_palm == ind_lup_pebblest) THEN |
| 3629 | luu_dep = 9 |
| 3630 | ELSEIF (luu_palm == ind_lup_woodchips) THEN |
| 3631 | luu_dep = 9 |
| 3632 | ELSEIF (luu_palm == ind_lup_tartan) THEN |
| 3633 | luu_dep = 9 |
| 3634 | ELSEIF (luu_palm == ind_lup_art_turf) THEN |
| 3635 | luu_dep = 9 |
| 3636 | ELSEIF (luu_palm == ind_lup_clay) THEN |
| 3637 | luu_dep = 9 |
| 3638 | ENDIF |
| 3639 | ENDIF |
| 3640 | |
| 3641 | IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN |
| 3642 | ! For windows in USM assume metal as this is |
| 3643 | ! as close as we get, assumed LU class desert |
| 3644 | ! as also assumed for pavements in LSM |
| 3645 | lud_palm = ind_lup_metal |
| 3646 | IF (lud_palm == ind_lup_user) THEN |
| 3647 | message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation' |
| 3648 | CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 ) |
| 3649 | ELSEIF (lud_palm == ind_lup_asph_conc) THEN |
| 3650 | lud_dep = 9 |
| 3651 | ELSEIF (lud_palm == ind_lup_asph) THEN |
| 3652 | lud_dep = 9 |
| 3653 | ELSEIF (lud_palm == ind_lup_conc) THEN |
| 3654 | lud_dep = 9 |
| 3655 | ELSEIF (lud_palm == ind_lup_sett) THEN |
| 3656 | lud_dep = 9 |
| 3657 | ELSEIF (lud_palm == ind_lup_pav_stones) THEN |
| 3658 | lud_dep = 9 |
| 3659 | ELSEIF (lud_palm == ind_lup_cobblest) THEN |
| 3660 | lud_dep = 9 |
| 3661 | ELSEIF (lud_palm == ind_lup_metal) THEN |
| 3662 | lud_dep = 9 |
| 3663 | ELSEIF (lud_palm == ind_lup_wood) THEN |
| 3664 | lud_dep = 9 |
| 3665 | ELSEIF (lud_palm == ind_lup_gravel) THEN |
| 3666 | lud_dep = 9 |
| 3667 | ELSEIF (lud_palm == ind_lup_f_gravel) THEN |
| 3668 | lud_dep = 9 |
| 3669 | ELSEIF (lud_palm == ind_lup_pebblest) THEN |
| 3670 | lud_dep = 9 |
| 3671 | ELSEIF (lud_palm == ind_lup_woodchips) THEN |
| 3672 | lud_dep = 9 |
| 3673 | ELSEIF (lud_palm == ind_lup_tartan) THEN |
| 3674 | lud_dep = 9 |
| 3675 | ELSEIF (lud_palm == ind_lup_art_turf) THEN |
| 3676 | lud_dep = 9 |
| 3677 | ELSEIF (lud_palm == ind_lup_clay) THEN |
| 3678 | lud_dep = 9 |
| 3679 | ENDIF |
| 3680 | ENDIF |
| 3681 | |
| 3682 | |
| 3683 | ! TODO: Activate these lines as soon as new ebsolver branch is merged: |
| 3684 | ! Set wetness indicator to dry or wet for usm vegetation or pavement |
| 3685 | !IF (surf_usm_h%c_liq(m) > 0) THEN |
| 3686 | ! nwet = 1 |
| 3687 | !ELSE |
| 3688 | nwet = 0 |
| 3689 | !ENDIF |
| 3690 | |
| 3691 | ! Compute length of time step |
| 3692 | IF ( call_chem_at_all_substeps ) THEN |
| 3693 | dt_chem = dt_3d * weight_pres(intermediate_timestep_count) |
| 3694 | ELSE |
| 3695 | dt_chem = dt_3d |
| 3696 | ENDIF |
| 3697 | |
| 3698 | |
| 3699 | dh = dzw(k) |
| 3700 | inv_dh = 1.0_wp / dh |
| 3701 | dt_dh = dt_chem/dh |
| 3702 | |
| 3703 | ! Concentration at i,j,k |
| 3704 | DO lsp = 1, nspec |
| 3705 | cc(lsp) = chem_species(lsp)%conc(k,j,i) |
| 3706 | ENDDO |
| 3707 | |
| 3708 | ! Temperature at i,j,k |
| 3709 | ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp |
| 3710 | |
| 3711 | ! Surface temperature in degrees Celcius: |
| 3712 | ts = ttemp - 273.15 ! in degrees celcius |
| 3713 | |
| 3714 | ! Viscosity of air in lowest layer |
| 3715 | visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) |
| 3716 | |
| 3717 | ! Air density in lowest layer |
| 3718 | dens = rho_air_zw(k) |
| 3719 | |
| 3720 | ! Calculate relative humidity from specific humidity for DEPAC |
| 3721 | qv_tmp = max(q(k,j,i),0.0_wp) |
| 3722 | relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) ) |
| 3723 | |
| 3724 | |
| 3725 | |
| 3726 | ! Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget |
| 3727 | ! for each surface fraction. Then derive overall budget taking into account the surface fractions. |
| 3728 | |
| 3729 | ! Walls/roofs |
| 3730 | IF (surf_usm_h%frac(ind_veg_wall,m) > 0) THEN |
| 3731 | |
| 3732 | |
| 3733 | ! No vegetation on non-green walls: |
| 3734 | lai = 0.0_wp |
| 3735 | sai = 0.0_wp |
| 3736 | |
| 3737 | slinnfac = 1.0_wp |
| 3738 | |
| 3739 | ! Get vd |
| 3740 | DO lsp = 1, nvar |
| 3741 | !Initialize |
| 3742 | vs = 0.0_wp |
| 3743 | vd_lu = 0.0_wp |
| 3744 | Rs = 0.0_wp |
| 3745 | Rb = 0.0_wp |
| 3746 | Rc_tot = 0.0_wp |
| 3747 | IF (spc_names(lsp) == 'PM10') THEN |
| 3748 | part_type = 1 |
| 3749 | ! sedimentation velocity in lowest layer |
| 3750 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3751 | particle_pars(ind_p_size, part_type), & |
| 3752 | particle_pars(ind_p_slip, part_type), & |
| 3753 | visc) |
| 3754 | |
| 3755 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3756 | vs, & |
| 3757 | particle_pars(ind_p_size, part_type), & |
| 3758 | particle_pars(ind_p_slip, part_type), & |
| 3759 | nwet, ttemp, dens, visc, & |
| 3760 | luu_dep, & |
| 3761 | r_aero_lu, ustar_lu) |
| 3762 | |
| 3763 | bud_now_luu(lsp) = - cc(lsp) * & |
| 3764 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3765 | |
| 3766 | |
| 3767 | ELSEIF (spc_names(lsp) == 'PM25') THEN |
| 3768 | part_type = 2 |
| 3769 | ! sedimentation velocity in lowest layer: |
| 3770 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3771 | particle_pars(ind_p_size, part_type), & |
| 3772 | particle_pars(ind_p_slip, part_type), & |
| 3773 | visc) |
| 3774 | |
| 3775 | |
| 3776 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3777 | vs, & |
| 3778 | particle_pars(ind_p_size, part_type), & |
| 3779 | particle_pars(ind_p_slip, part_type), & |
| 3780 | nwet, ttemp, dens, visc, & |
| 3781 | luu_dep , & |
| 3782 | r_aero_lu, ustar_lu) |
| 3783 | |
| 3784 | bud_now_luu(lsp) = - cc(lsp) * & |
| 3785 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3786 | |
| 3787 | |
| 3788 | ELSE |
| 3789 | |
| 3790 | ! GASES |
| 3791 | |
| 3792 | ! Read spc_name of current species for gas parameter |
| 3793 | |
| 3794 | IF (ANY(pspecnames(:) == spc_names(lsp))) THEN |
| 3795 | i_pspec = 0 |
| 3796 | DO pspec = 1, 69 |
| 3797 | IF (pspecnames(pspec) == spc_names(lsp)) THEN |
| 3798 | i_pspec = pspec |
| 3799 | END IF |
| 3800 | ENDDO |
| 3801 | |
| 3802 | ELSE |
| 3803 | ! Species for now not deposited |
| 3804 | CYCLE |
| 3805 | ENDIF |
| 3806 | |
| 3807 | ! Factor used for conversion from ppb to ug/m3 : |
| 3808 | ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & |
| 3809 | ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) |
| 3810 | ! c 1e-9 xm_tracer 1e9 / xm_air dens |
| 3811 | ! thus: |
| 3812 | ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 |
| 3813 | ! Use density at lowest layer: |
| 3814 | |
| 3815 | ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 |
| 3816 | |
| 3817 | ! atmospheric concentration in DEPAC is requested in ug/m3: |
| 3818 | ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole |
| 3819 | catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 |
| 3820 | |
| 3821 | ! Diffusivity for DEPAC relevant gases |
| 3822 | ! Use default value |
| 3823 | diffc = 0.11e-4 |
| 3824 | ! overwrite with known coefficients of diffusivity from Massman (1998) |
| 3825 | IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 |
| 3826 | IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 |
| 3827 | IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 |
| 3828 | IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 |
| 3829 | IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 |
| 3830 | IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 |
| 3831 | IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 |
| 3832 | |
| 3833 | |
| 3834 | ! Get quasi-laminar boundary layer resistance Rb: |
| 3835 | CALL get_rb_cell( (luu_dep == ilu_water_sea) .or. (luu_dep == ilu_water_inland), & |
| 3836 | z0h_lu, ustar_lu, diffc, & |
| 3837 | Rb) |
| 3838 | |
| 3839 | ! Get Rc_tot |
| 3840 | CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & |
| 3841 | relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & |
| 3842 | r_aero_lu , Rb) |
| 3843 | |
| 3844 | |
| 3845 | ! Calculate budget |
| 3846 | IF (Rc_tot .le. 0.0) THEN |
| 3847 | |
| 3848 | bud_now_luu(lsp) = 0.0_wp |
| 3849 | |
| 3850 | ELSE |
| 3851 | |
| 3852 | ! Compute exchange velocity for current lu: |
| 3853 | vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) |
| 3854 | |
| 3855 | bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & |
| 3856 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3857 | ENDIF |
| 3858 | |
| 3859 | ENDIF |
| 3860 | ENDDO |
| 3861 | ENDIF |
| 3862 | |
| 3863 | |
| 3864 | ! Green usm surfaces |
| 3865 | IF (surf_usm_h%frac(ind_pav_green,m) > 0) THEN |
| 3866 | |
| 3867 | |
| 3868 | slinnfac = 1.0_wp |
| 3869 | |
| 3870 | ! Get vd |
| 3871 | DO lsp = 1, nvar |
| 3872 | !Initialize |
| 3873 | vs = 0.0_wp |
| 3874 | vd_lu = 0.0_wp |
| 3875 | Rs = 0.0_wp |
| 3876 | Rb = 0.0_wp |
| 3877 | Rc_tot = 0.0_wp |
| 3878 | IF (spc_names(lsp) == 'PM10') THEN |
| 3879 | part_type = 1 |
| 3880 | ! sedimentation velocity in lowest layer: |
| 3881 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3882 | particle_pars(ind_p_size, part_type), & |
| 3883 | particle_pars(ind_p_slip, part_type), & |
| 3884 | visc) |
| 3885 | |
| 3886 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3887 | vs, & |
| 3888 | particle_pars(ind_p_size, part_type), & |
| 3889 | particle_pars(ind_p_slip, part_type), & |
| 3890 | nwet, ttemp, dens, visc, & |
| 3891 | lug_dep, & |
| 3892 | r_aero_lu, ustar_lu) |
| 3893 | |
| 3894 | bud_now_lug(lsp) = - cc(lsp) * & |
| 3895 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3896 | |
| 3897 | |
| 3898 | ELSEIF (spc_names(lsp) == 'PM25') THEN |
| 3899 | part_type = 2 |
| 3900 | ! sedimentation velocity in lowest layer: |
| 3901 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 3902 | particle_pars(ind_p_size, part_type), & |
| 3903 | particle_pars(ind_p_slip, part_type), & |
| 3904 | visc) |
| 3905 | |
| 3906 | |
| 3907 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 3908 | vs, & |
| 3909 | particle_pars(ind_p_size, part_type), & |
| 3910 | particle_pars(ind_p_slip, part_type), & |
| 3911 | nwet, ttemp, dens, visc, & |
| 3912 | lug_dep, & |
| 3913 | r_aero_lu, ustar_lu) |
| 3914 | |
| 3915 | bud_now_lug(lsp) = - cc(lsp) * & |
| 3916 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3917 | |
| 3918 | |
| 3919 | ELSE |
| 3920 | |
| 3921 | ! GASES |
| 3922 | |
| 3923 | ! Read spc_name of current species for gas parameter |
| 3924 | |
| 3925 | IF (ANY(pspecnames(:) == spc_names(lsp))) THEN |
| 3926 | i_pspec = 0 |
| 3927 | DO pspec = 1, 69 |
| 3928 | IF (pspecnames(pspec) == spc_names(lsp)) THEN |
| 3929 | i_pspec = pspec |
| 3930 | END IF |
| 3931 | ENDDO |
| 3932 | |
| 3933 | ELSE |
| 3934 | ! Species for now not deposited |
| 3935 | CYCLE |
| 3936 | ENDIF |
| 3937 | |
| 3938 | ! Factor used for conversion from ppb to ug/m3 : |
| 3939 | ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & |
| 3940 | ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) |
| 3941 | ! c 1e-9 xm_tracer 1e9 / xm_air dens |
| 3942 | ! thus: |
| 3943 | ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 |
| 3944 | ! Use density at lowest layer: |
| 3945 | |
| 3946 | ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 |
| 3947 | |
| 3948 | ! atmospheric concentration in DEPAC is requested in ug/m3: |
| 3949 | ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole |
| 3950 | catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 |
| 3951 | |
| 3952 | ! Diffusivity for DEPAC relevant gases |
| 3953 | ! Use default value |
| 3954 | diffc = 0.11e-4 |
| 3955 | ! overwrite with known coefficients of diffusivity from Massman (1998) |
| 3956 | IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 |
| 3957 | IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 |
| 3958 | IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 |
| 3959 | IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 |
| 3960 | IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 |
| 3961 | IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 |
| 3962 | IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 |
| 3963 | |
| 3964 | |
| 3965 | ! Get quasi-laminar boundary layer resistance Rb: |
| 3966 | CALL get_rb_cell( (lug_dep == ilu_water_sea) .or. (lug_dep == ilu_water_inland), & |
| 3967 | z0h_lu, ustar_lu, diffc, & |
| 3968 | Rb) |
| 3969 | |
| 3970 | |
| 3971 | ! Get Rc_tot |
| 3972 | CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & |
| 3973 | relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & |
| 3974 | r_aero_lu , Rb) |
| 3975 | |
| 3976 | |
| 3977 | ! Calculate budget |
| 3978 | IF (Rc_tot .le. 0.0) THEN |
| 3979 | |
| 3980 | bud_now_lug(lsp) = 0.0_wp |
| 3981 | |
| 3982 | ELSE |
| 3983 | |
| 3984 | ! Compute exchange velocity for current lu: |
| 3985 | vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) |
| 3986 | |
| 3987 | bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & |
| 3988 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 3989 | ENDIF |
| 3990 | |
| 3991 | |
| 3992 | ENDIF |
| 3993 | ENDDO |
| 3994 | ENDIF |
| 3995 | |
| 3996 | |
| 3997 | ! Windows |
| 3998 | IF (surf_usm_h%frac(ind_wat_win,m) > 0) THEN |
| 3999 | |
| 4000 | |
| 4001 | ! No vegetation on windows: |
| 4002 | lai = 0.0_wp |
| 4003 | sai = 0.0_wp |
| 4004 | |
| 4005 | slinnfac = 1.0_wp |
| 4006 | |
| 4007 | ! Get vd |
| 4008 | DO lsp = 1, nvar |
| 4009 | !Initialize |
| 4010 | vs = 0.0_wp |
| 4011 | vd_lu = 0.0_wp |
| 4012 | Rs = 0.0_wp |
| 4013 | Rb = 0.0_wp |
| 4014 | Rc_tot = 0.0_wp |
| 4015 | IF (spc_names(lsp) == 'PM10') THEN |
| 4016 | part_type = 1 |
| 4017 | ! sedimentation velocity in lowest layer: |
| 4018 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 4019 | particle_pars(ind_p_size, part_type), & |
| 4020 | particle_pars(ind_p_slip, part_type), & |
| 4021 | visc) |
| 4022 | |
| 4023 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 4024 | vs, & |
| 4025 | particle_pars(ind_p_size, part_type), & |
| 4026 | particle_pars(ind_p_slip, part_type), & |
| 4027 | nwet, ttemp, dens, visc, & |
| 4028 | lud_dep, & |
| 4029 | r_aero_lu, ustar_lu) |
| 4030 | |
| 4031 | bud_now_lud(lsp) = - cc(lsp) * & |
| 4032 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 4033 | |
| 4034 | |
| 4035 | ELSEIF (spc_names(lsp) == 'PM25') THEN |
| 4036 | part_type = 2 |
| 4037 | ! sedimentation velocity in lowest layer: |
| 4038 | vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & |
| 4039 | particle_pars(ind_p_size, part_type), & |
| 4040 | particle_pars(ind_p_slip, part_type), & |
| 4041 | visc) |
| 4042 | |
| 4043 | |
| 4044 | CALL drydepo_aero_zhang_vd( vd_lu, Rs, & |
| 4045 | vs, & |
| 4046 | particle_pars(ind_p_size, part_type), & |
| 4047 | particle_pars(ind_p_slip, part_type), & |
| 4048 | nwet, ttemp, dens, visc, & |
| 4049 | lud_dep, & |
| 4050 | r_aero_lu, ustar_lu) |
| 4051 | |
| 4052 | bud_now_lud(lsp) = - cc(lsp) * & |
| 4053 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 4054 | |
| 4055 | |
| 4056 | ELSE |
| 4057 | |
| 4058 | ! GASES |
| 4059 | |
| 4060 | ! Read spc_name of current species for gas PARAMETER |
| 4061 | |
| 4062 | IF (ANY(pspecnames(:) == spc_names(lsp))) THEN |
| 4063 | i_pspec = 0 |
| 4064 | DO pspec = 1, 69 |
| 4065 | IF (pspecnames(pspec) == spc_names(lsp)) THEN |
| 4066 | i_pspec = pspec |
| 4067 | END IF |
| 4068 | ENDDO |
| 4069 | |
| 4070 | ELSE |
| 4071 | ! Species for now not deposited |
| 4072 | CYCLE |
| 4073 | ENDIF |
| 4074 | |
| 4075 | ! Factor used for conversion from ppb to ug/m3 : |
| 4076 | ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & |
| 4077 | ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) |
| 4078 | ! c 1e-9 xm_tracer 1e9 / xm_air dens |
| 4079 | ! thus: |
| 4080 | ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 |
| 4081 | ! Use density at lowest layer: |
| 4082 | |
| 4083 | ppm_to_ugm3 = (dens/xm_air)*0.001_wp ! (mole air)/m3 |
| 4084 | |
| 4085 | ! atmospheric concentration in DEPAC is requested in ug/m3: |
| 4086 | ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole |
| 4087 | catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 |
| 4088 | |
| 4089 | ! Diffusivity for DEPAC relevant gases |
| 4090 | ! Use default value |
| 4091 | diffc = 0.11e-4 |
| 4092 | ! overwrite with known coefficients of diffusivity from Massman (1998) |
| 4093 | IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 |
| 4094 | IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 |
| 4095 | IF ( spc_names(lsp) == 'O3' ) diffc = 0.144e-4 |
| 4096 | IF ( spc_names(lsp) == 'CO' ) diffc = 0.176e-4 |
| 4097 | IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4 |
| 4098 | IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 |
| 4099 | IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 |
| 4100 | |
| 4101 | |
| 4102 | ! Get quasi-laminar boundary layer resistance Rb: |
| 4103 | CALL get_rb_cell( (lud_dep == ilu_water_sea) .or. (lud_dep == ilu_water_inland), & |
| 4104 | z0h_lu, ustar_lu, diffc, & |
| 4105 | Rb) |
| 4106 | |
| 4107 | ! Get Rc_tot |
| 4108 | CALL drydepos_gas_depac(spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), & |
| 4109 | relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & |
| 4110 | r_aero_lu , Rb) |
| 4111 | |
| 4112 | |
| 4113 | ! Calculate budget |
| 4114 | IF (Rc_tot .le. 0.0) THEN |
| 4115 | |
| 4116 | bud_now_lud(lsp) = 0.0_wp |
| 4117 | |
| 4118 | ELSE |
| 4119 | |
| 4120 | ! Compute exchange velocity for current lu: |
| 4121 | vd_lu = 1.0 / (r_aero_lu + Rb + Rc_tot ) |
| 4122 | |
| 4123 | bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & |
| 4124 | (1.0 - exp(-vd_lu*dt_dh ))*dh |
| 4125 | ENDIF |
| 4126 | |
| 4127 | ENDIF |
| 4128 | ENDDO |
| 4129 | ENDIF |
| 4130 | |
| 4131 | |
| 4132 | bud_now = 0.0_wp |
| 4133 | ! Calculate budget for surface m and adapt concentration |
| 4134 | DO lsp = 1, nspec |
| 4135 | |
| 4136 | |
| 4137 | bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m)*bud_now_luu(lsp) + & |
| 4138 | surf_usm_h%frac(ind_pav_green,m)*bud_now_lug(lsp) + & |
| 4139 | surf_usm_h%frac(ind_wat_win,m)*bud_now_lud(lsp) |
| 4140 | |
| 4141 | ! Compute new concentration, add contribution of all exchange processes: |
| 4142 | cc(lsp) = cc(lsp) + bud_now(lsp)*inv_dh |
| 4143 | |
| 4144 | chem_species(lsp)%conc(k,j,i) = max(0.0_wp, cc(lsp)) |
| 4145 | |
| 4146 | ENDDO |
| 4147 | |
| 4148 | ENDIF |
| 4149 | |
| 4150 | |
| 4151 | |
| 4152 | END SUBROUTINE chem_depo |
| 4153 | |
| 4154 | |
| 4155 | |
| 4156 | !---------------------------------------------------------------------------------- |
| 4157 | ! |
| 4158 | ! DEPAC: |
| 4159 | ! Code of the DEPAC routine and corresponding subroutines below from the DEPAC |
| 4160 | ! module of the LOTOS-EUROS model (Manders et al., 2017) |
| 4161 | ! |
| 4162 | ! Original DEPAC routines by RIVM and TNO (2015), for Documentation see |
| 4163 | ! van Zanten et al., 2010. |
| 4164 | !--------------------------------------------------------------------------------- |
| 4165 | ! |
| 4166 | !---------------------------------------------------------------------------------- |
| 4167 | ! |
| 4168 | ! depac: compute total canopy (or surface) resistance Rc for gases |
| 4169 | !---------------------------------------------------------------------------------- |
| 4170 | |
| 4171 | SUBROUTINE drydepos_gas_depac(compnam, day_of_year, lat, t, ust, glrad, sinphi, & |
| 4172 | rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc, & |
| 4173 | ra, rb) |
| 4174 | |
| 4175 | |
| 4176 | ! The last four rows of depac arguments are OPTIONAL: |
| 4177 | ! |
| 4178 | ! A. compute Rc_tot without compensation points (ccomp_tot will be zero): |
| 4179 | ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi]) |
| 4180 | ! |
| 4181 | ! B. compute Rc_tot with compensation points (used for LOTOS-EUROS): |
| 4182 | ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & |
| 4183 | ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water) |
| 4184 | ! |
| 4185 | ! C. compute effective Rc based on compensation points (used for OPS): |
| 4186 | ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & |
| 4187 | ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & |
| 4188 | ! ra, rb, rc_eff) |
| 4189 | ! X1. Extra (OPTIONAL) output variables: |
| 4190 | ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & |
| 4191 | ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & |
| 4192 | ! ra, rb, rc_eff, & |
| 4193 | ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out) |
| 4194 | ! X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): |
| 4195 | ! CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & |
| 4196 | ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & |
| 4197 | ! ra, rb, rc_eff, & |
| 4198 | ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, & |
| 4199 | ! calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA) |
| 4200 | |
| 4201 | |
| 4202 | IMPLICIT NONE |
| 4203 | |
| 4204 | CHARACTER(len=*) , INTENT(in) :: compnam ! component name |
| 4205 | ! 'HNO3','NO','NO2','O3','SO2','NH3' |
| 4206 | INTEGER(iwp) , INTENT(in) :: day_of_year ! day of year, 1 ... 365 (366) |
| 4207 | REAL(kind=wp) , INTENT(in) :: lat ! latitude Northern hemisphere (degrees) (DEPAC cannot be used for S. hemisphere) |
| 4208 | REAL(kind=wp) , INTENT(in) :: t ! temperature (C) |
| 4209 | ! NB discussion issue is temp T_2m or T_surf or T_leaf? |
| 4210 | REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) |
| 4211 | REAL(kind=wp) , INTENT(in) :: glrad ! global radiation (W/m2) |
| 4212 | REAL(kind=wp) , INTENT(in) :: sinphi ! sin of solar elevation angle |
| 4213 | REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) |
| 4214 | REAL(kind=wp) , INTENT(in) :: lai ! one-sidedleaf area index (-) |
| 4215 | REAL(kind=wp) , INTENT(in) :: sai ! surface area index (-) (lai + branches and stems) |
| 4216 | REAL(kind=wp) , INTENT(in) :: diffc ! Diffusivity |
| 4217 | INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow |
| 4218 | INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu |
| 4219 | INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio used for SO2; |
| 4220 | ! iratns = 1: low NH3/SO2 |
| 4221 | ! iratns = 2: high NH3/SO2 |
| 4222 | ! iratns = 3: very low NH3/SO2 |
| 4223 | REAL(kind=wp) , INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) |
| 4224 | REAL(kind=wp) , INTENT(out) :: ccomp_tot ! total compensation point (ug/m3) [= 0 for species that don't have a compensation |
| 4225 | REAL(kind=wp) ,INTENT(in) :: p ! pressure (Pa) |
| 4226 | REAL(kind=wp), INTENT(in) :: catm ! actual atmospheric concentration (ug/m3) |
| 4227 | REAL(kind=wp), INTENT(in) :: ra ! aerodynamic resistance (s/m) |
| 4228 | REAL(kind=wp), INTENT(in) :: rb ! boundary layer resistance (s/m) |
| 4229 | |
| 4230 | ! Local variables: |
| 4231 | REAL(kind=wp) :: laimax ! maximum leaf area index (-) |
| 4232 | LOGICAL :: ready ! Rc has been set |
| 4233 | ! = 1 -> constant Rc |
| 4234 | ! = 2 -> temperature dependent Rc |
| 4235 | REAL(kind=wp) :: gw ! external leaf conductance (m/s) |
| 4236 | REAL(kind=wp) :: gstom ! stomatal conductance (m/s) |
| 4237 | REAL(kind=wp) :: gsoil_eff ! effective soil conductance (m/s) |
| 4238 | REAL(kind=wp) :: gc_tot ! total canopy conductance (m/s) |
| 4239 | REAL(kind=wp) :: cw ! external leaf surface compensation point (ug/m3) |
| 4240 | REAL(kind=wp) :: cstom ! stomatal compensation point (ug/m3) |
| 4241 | REAL(kind=wp) :: csoil ! soil compensation point (ug/m3) |
| 4242 | |
| 4243 | ! Vegetation indicators: |
| 4244 | LOGICAL :: LAI_present ! leaves are present for current land use type |
| 4245 | LOGICAL :: SAI_present ! vegetation is present for current land use type |
| 4246 | |
| 4247 | ! Component number taken from component name, paramteres matched with include files |
| 4248 | INTEGER(iwp) :: icmp |
| 4249 | |
| 4250 | ! component numbers: |
| 4251 | INTEGER(iwp), PARAMETER :: icmp_o3 = 1 |
| 4252 | INTEGER(iwp), PARAMETER :: icmp_so2 = 2 |
| 4253 | INTEGER(iwp), PARAMETER :: icmp_no2 = 3 |
| 4254 | INTEGER(iwp), PARAMETER :: icmp_no = 4 |
| 4255 | INTEGER(iwp), PARAMETER :: icmp_nh3 = 5 |
| 4256 | INTEGER(iwp), PARAMETER :: icmp_co = 6 |
| 4257 | INTEGER(iwp), PARAMETER :: icmp_no3 = 7 |
| 4258 | INTEGER(iwp), PARAMETER :: icmp_hno3 = 8 |
| 4259 | INTEGER(iwp), PARAMETER :: icmp_n2o5 = 9 |
| 4260 | INTEGER(iwp), PARAMETER :: icmp_h2o2 = 10 |
| 4261 | |
| 4262 | |
| 4263 | |
| 4264 | ! Define component number |
| 4265 | SELECT CASE ( trim(compnam) ) |
| 4266 | |
| 4267 | CASE ( 'O3', 'o3' ) |
| 4268 | icmp = icmp_o3 |
| 4269 | |
| 4270 | CASE ( 'SO2', 'so2' ) |
| 4271 | icmp = icmp_so2 |
| 4272 | |
| 4273 | CASE ( 'NO2', 'no2' ) |
| 4274 | icmp = icmp_no2 |
| 4275 | |
| 4276 | CASE ( 'NO', 'no' ) |
| 4277 | icmp = icmp_no |
| 4278 | |
| 4279 | CASE ( 'NH3', 'nh3' ) |
| 4280 | icmp = icmp_nh3 |
| 4281 | |
| 4282 | CASE ( 'CO', 'co' ) |
| 4283 | icmp = icmp_co |
| 4284 | |
| 4285 | CASE ( 'NO3', 'no3' ) |
| 4286 | icmp = icmp_no3 |
| 4287 | |
| 4288 | CASE ( 'HNO3', 'hno3' ) |
| 4289 | icmp = icmp_hno3 |
| 4290 | |
| 4291 | CASE ( 'N2O5', 'n2o5' ) |
| 4292 | icmp = icmp_n2o5 |
| 4293 | |
| 4294 | CASE ( 'H2O2', 'h2o2' ) |
| 4295 | icmp = icmp_h2o2 |
| 4296 | |
| 4297 | CASE default |
| 4298 | !Component not part of DEPAC --> not deposited |
| 4299 | RETURN |
| 4300 | |
| 4301 | END SELECT |
| 4302 | |
| 4303 | ! inititalize |
| 4304 | gw = 0. |
| 4305 | gstom = 0. |
| 4306 | gsoil_eff = 0. |
| 4307 | gc_tot = 0. |
| 4308 | cw = 0. |
| 4309 | cstom = 0. |
| 4310 | csoil = 0. |
| 4311 | |
| 4312 | |
| 4313 | ! Check whether vegetation is present (in that CASE the leaf or surface area index > 0): |
| 4314 | LAI_present = (lai .gt. 0.0) |
| 4315 | SAI_present = (sai .gt. 0.0) |
| 4316 | |
| 4317 | ! Set Rc (i.e. rc_tot) in special cases: |
| 4318 | CALL rc_special(icmp,compnam,lu,t, nwet,rc_tot,ready,ccomp_tot) |
| 4319 | |
| 4320 | ! set effective resistance equal to total resistance, this will be changed for compensation points |
| 4321 | !IF ( present(rc_eff) ) then |
| 4322 | ! rc_eff = rc_tot |
| 4323 | !END IF |
| 4324 | |
| 4325 | ! IF Rc is not set: |
| 4326 | IF (.not. ready) then |
| 4327 | |
| 4328 | ! External conductance: |
| 4329 | CALL rc_gw(compnam,iratns,t,rh,nwet,SAI_present,sai,gw) |
| 4330 | |
| 4331 | ! Stomatal conductance: |
| 4332 | ! CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p,& |
| 4333 | ! smi=smi,calc_stom_o3flux=calc_stom_o3flux) |
| 4334 | CALL rc_gstom(icmp,compnam,lu,LAI_present,lai,glrad,sinphi,t,rh,diffc,gstom,p) |
| 4335 | ! Effective soil conductance: |
| 4336 | CALL rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff) |
| 4337 | |
| 4338 | ! Total canopy conductance (gc_tot) and resistance Rc (rc_tot): |
| 4339 | CALL rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot) |
| 4340 | |
| 4341 | ! Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3): |
| 4342 | ! CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,& |
| 4343 | ! LAI_present, SAI_present, & |
| 4344 | ! ccomp_tot, & |
| 4345 | ! catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, & |
| 4346 | ! c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, & |
| 4347 | ! tsea=tsea,cw=cw,cstom=cstom,csoil=csoil ) |
| 4348 | |
| 4349 | ! Effective Rc based on compensation points: |
| 4350 | ! IF ( present(rc_eff) ) then |
| 4351 | ! check on required arguments: |
| 4352 | !IF ( (.not. present(catm)) .or. (.not. present(ra)) .or. (.not. present(rb)) ) then |
| 4353 | ! stop 'output argument rc_eff requires input arguments catm, ra and rb' |
| 4354 | !END IF |
| 4355 | ! compute rc_eff : |
| 4356 | !CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) |
| 4357 | !ENDIF |
| 4358 | ENDIF |
| 4359 | |
| 4360 | END SUBROUTINE drydepos_gas_depac |
| 4361 | |
| 4362 | |
| 4363 | |
| 4364 | !------------------------------------------------------------------- |
| 4365 | ! rc_special: compute total canopy resistance in special CASEs |
| 4366 | !------------------------------------------------------------------- |
| 4367 | SUBROUTINE rc_special(icmp,compnam,lu,t,nwet,rc_tot,ready,ccomp_tot) |
| 4368 | |
| 4369 | INTEGER(iwp) , INTENT(in) :: icmp ! component index |
| 4370 | CHARACTER(len=*) , INTENT(in) :: compnam ! component name |
| 4371 | INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu |
| 4372 | REAL(kind=wp) , INTENT(in) :: t ! temperature (C) |
| 4373 | ! = 1 -> constant Rc |
| 4374 | ! = 2 -> temperature dependent Rc |
| 4375 | INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow |
| 4376 | REAL(kind=wp) , INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) |
| 4377 | LOGICAL , INTENT(out) :: ready ! Rc has been set |
| 4378 | REAL(kind=wp) , INTENT(out) :: ccomp_tot ! total compensation point (ug/m3) |
| 4379 | |
| 4380 | ! rc_tot is not yet set: |
| 4381 | ready = .false. |
| 4382 | |
| 4383 | ! Default compensation point in special CASEs = 0: |
| 4384 | ccomp_tot = 0.0 |
| 4385 | |
| 4386 | SELECT CASE(trim(compnam)) |
| 4387 | CASE('HNO3','N2O5','NO3','H2O2') |
| 4388 | ! No separate resistances for HNO3; just one total canopy resistance: |
| 4389 | IF (t .lt. -5.0 .and. nwet .eq. 9) then |
| 4390 | ! T < 5 C and snow: |
| 4391 | rc_tot = 50. |
| 4392 | ELSE |
| 4393 | ! all other circumstances: |
| 4394 | rc_tot = 10.0 |
| 4395 | ENDIF |
| 4396 | ready = .true. |
| 4397 | |
| 4398 | CASE('NO','CO') |
| 4399 | IF (lu .eq. ilu_water_sea .or. lu .eq. ilu_water_inland) then ! water |
| 4400 | rc_tot = 2000. |
| 4401 | ready = .true. |
| 4402 | ELSEIF (nwet .eq. 1) then ! wet |
| 4403 | rc_tot = 2000. |
| 4404 | ready = .true. |
| 4405 | ENDIF |
| 4406 | CASE('NO2','O3','SO2','NH3') |
| 4407 | ! snow surface: |
| 4408 | IF (nwet.eq.9) then |
| 4409 | !CALL rc_snow(ipar_snow(icmp),t,rc_tot) |
| 4410 | ready = .true. |
| 4411 | ENDIF |
| 4412 | CASE default |
| 4413 | message_string = 'Component '// trim(compnam) // ' not supported' |
| 4414 | CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 ) |
| 4415 | END SELECT |
| 4416 | |
| 4417 | END SUBROUTINE rc_special |
| 4418 | |
| 4419 | |
| 4420 | |
| 4421 | !------------------------------------------------------------------- |
| 4422 | ! rc_gw: compute external conductance |
| 4423 | !------------------------------------------------------------------- |
| 4424 | SUBROUTINE rc_gw( compnam, iratns,t,rh,nwet, SAI_present, sai, gw ) |
| 4425 | |
| 4426 | ! Input/output variables: |
| 4427 | CHARACTER(len=*) , INTENT(in) :: compnam ! component name |
| 4428 | INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio; |
| 4429 | ! iratns = 1: low NH3/SO2 |
| 4430 | ! iratns = 2: high NH3/SO2 |
| 4431 | ! iratns = 3: very low NH3/SO2 |
| 4432 | REAL(kind=wp) , INTENT(in) :: t ! temperature (C) |
| 4433 | REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) |
| 4434 | INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow |
| 4435 | LOGICAL , INTENT(in) :: SAI_present |
| 4436 | REAL(kind=wp) , INTENT(in) :: sai ! one-sided leaf area index (-) |
| 4437 | REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) |
| 4438 | |
| 4439 | SELECT CASE(trim(compnam)) |
| 4440 | |
| 4441 | !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3 |
| 4442 | |
| 4443 | CASE('NO2') |
| 4444 | CALL rw_constant(2000.0_wp,SAI_present,gw) |
| 4445 | |
| 4446 | CASE('NO','CO') |
| 4447 | CALL rw_constant(-9999.0_wp,SAI_present,gw) ! see Erisman et al, 1994 section 3.2.3 |
| 4448 | |
| 4449 | CASE('O3') |
| 4450 | CALL rw_constant(2500.0_wp,SAI_present,gw) |
| 4451 | |
| 4452 | CASE('SO2') |
| 4453 | CALL rw_so2( t, nwet, rh, iratns, SAI_present, gw ) |
| 4454 | |
| 4455 | CASE('NH3') |
| 4456 | CALL rw_nh3_sutton(t,rh,SAI_present,gw) |
| 4457 | |
| 4458 | ! conversion from leaf resistance to canopy resistance by multiplying with SAI: |
| 4459 | Gw = sai*gw |
| 4460 | |
| 4461 | CASE default |
| 4462 | message_string = 'Component '// trim(compnam) // ' not supported' |
| 4463 | CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 ) |
| 4464 | END SELECT |
| 4465 | |
| 4466 | END SUBROUTINE rc_gw |
| 4467 | |
| 4468 | |
| 4469 | |
| 4470 | !------------------------------------------------------------------- |
| 4471 | ! rw_so2: compute external leaf conductance for SO2 |
| 4472 | !------------------------------------------------------------------- |
| 4473 | SUBROUTINE rw_so2( t, nwet, rh, iratns, SAI_present, gw ) |
| 4474 | |
| 4475 | ! Input/output variables: |
| 4476 | REAL(kind=wp) , INTENT(in) :: t ! temperature (C) |
| 4477 | INTEGER(iwp) , INTENT(in) :: nwet ! wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow |
| 4478 | REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) |
| 4479 | INTEGER(iwp) , INTENT(in) :: iratns ! index for NH3/SO2 ratio; |
| 4480 | ! iratns = 1: low NH3/SO2 |
| 4481 | ! iratns = 2: high NH3/SO2 |
| 4482 | ! iratns = 3: very low NH3/SO2 |
| 4483 | LOGICAL, INTENT(in) :: SAI_present |
| 4484 | REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) |
| 4485 | |
| 4486 | ! Variables from module: |
| 4487 | ! SAI_present: vegetation is present |
| 4488 | |
| 4489 | ! Local variables: |
| 4490 | REAL(kind=wp) :: rw ! external leaf resistance (s/m) |
| 4491 | |
| 4492 | ! Check IF vegetation present: |
| 4493 | IF (SAI_present) then |
| 4494 | |
| 4495 | IF (nwet .eq. 0) then |
| 4496 | !-------------------------- |
| 4497 | ! dry surface |
| 4498 | !-------------------------- |
| 4499 | ! T > -1 C |
| 4500 | IF (t .gt. -1.0) then |
| 4501 | IF (rh .lt. 81.3) then |
| 4502 | rw = 25000*exp(-0.0693*rh) |
| 4503 | ELSE |
| 4504 | rw = 0.58e12*exp(-0.278*rh) + 10. |
| 4505 | ENDIF |
| 4506 | ELSE |
| 4507 | ! -5 C < T <= -1 C |
| 4508 | IF (t .gt. -5.0) then |
| 4509 | rw=200 |
| 4510 | ELSE |
| 4511 | ! T <= -5 C |
| 4512 | rw=500 |
| 4513 | ENDIF |
| 4514 | ENDIF |
| 4515 | ELSE |
| 4516 | !-------------------------- |
| 4517 | ! wet surface |
| 4518 | !-------------------------- |
| 4519 | rw = 10. !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10 |
| 4520 | ENDIF |
| 4521 | |
| 4522 | ! very low NH3/SO2 ratio: |
| 4523 | IF (iratns == iratns_very_low ) rw = rw+50. |
| 4524 | |
| 4525 | ! Conductance: |
| 4526 | gw = 1./rw |
| 4527 | ELSE |
| 4528 | ! no vegetation: |
| 4529 | gw = 0.0 |
| 4530 | ENDIF |
| 4531 | |
| 4532 | END SUBROUTINE rw_so2 |
| 4533 | |
| 4534 | |
| 4535 | |
| 4536 | !------------------------------------------------------------------- |
| 4537 | ! rw_nh3_sutton: compute external leaf conductance for NH3, |
| 4538 | ! following Sutton & Fowler, 1993 |
| 4539 | !------------------------------------------------------------------- |
| 4540 | SUBROUTINE rw_nh3_sutton(tsurf,rh,SAI_present,gw) |
| 4541 | |
| 4542 | ! Input/output variables: |
| 4543 | REAL(kind=wp) , INTENT(in) :: tsurf ! surface temperature (C) |
| 4544 | REAL(kind=wp) , INTENT(in) :: rh ! relative humidity (%) |
| 4545 | LOGICAL, INTENT(in) :: SAI_present |
| 4546 | REAL(kind=wp) , INTENT(out) :: gw ! external leaf conductance (m/s) |
| 4547 | |
| 4548 | ! Variables from module: |
| 4549 | ! SAI_present: vegetation is present |
| 4550 | |
| 4551 | ! Local variables: |
| 4552 | REAL(kind=wp) :: rw ! external leaf resistance (s/m) |
| 4553 | REAL(kind=wp) :: sai_grass_haarweg ! surface area index at experimental site Haarweg |
| 4554 | |
| 4555 | ! Fix SAI_grass at value valid for Haarweg data for which gamma_w parametrization is derived |
| 4556 | sai_grass_haarweg = 3.5 |
| 4557 | |
| 4558 | ! 100 - rh |
| 4559 | ! rw = 2.0 * exp(----------) |
| 4560 | ! 12 |
| 4561 | |
| 4562 | IF (SAI_present) then |
| 4563 | |
| 4564 | ! External resistance according to Sutton & Fowler, 1993 |
| 4565 | rw = 2.0 * exp((100.0 - rh)/12.0) |
| 4566 | rw = sai_grass_haarweg * rw |
| 4567 | |
| 4568 | ! Frozen soil (from Depac v1): |
| 4569 | IF (tsurf .lt. 0) rw = 200 |
| 4570 | |
| 4571 | ! Conductance: |
| 4572 | gw = 1./rw |
| 4573 | ELSE |
| 4574 | ! no vegetation: |
| 4575 | gw = 0.0 |
| 4576 | ENDIF |
| 4577 | |
| 4578 | END SUBROUTINE rw_nh3_sutton |
| 4579 | |
| 4580 | |
| 4581 | |
| 4582 | !------------------------------------------------------------------- |
| 4583 | ! rw_constant: compute constant external leaf conductance |
| 4584 | !------------------------------------------------------------------- |
| 4585 | SUBROUTINE rw_constant(rw_val,SAI_present,gw) |
| 4586 | |
| 4587 | ! Input/output variables: |
| 4588 | REAL(kind=wp) , INTENT(in) :: rw_val ! constant value of Rw |
| 4589 | LOGICAL , INTENT(in) :: SAI_present |
| 4590 | REAL(kind=wp) , INTENT(out) :: gw ! wernal leaf conductance (m/s) |
| 4591 | |
| 4592 | ! Variables from module: |
| 4593 | ! SAI_present: vegetation is present |
| 4594 | |
| 4595 | ! Compute conductance: |
| 4596 | IF (SAI_present .and. .not.missing(rw_val)) then |
| 4597 | gw = 1./rw_val |
| 4598 | ELSE |
| 4599 | gw = 0. |
| 4600 | ENDIF |
| 4601 | |
| 4602 | END SUBROUTINE rw_constant |
| 4603 | |
| 4604 | |
| 4605 | |
| 4606 | !------------------------------------------------------------------- |
| 4607 | ! rc_gstom: compute stomatal conductance |
| 4608 | !------------------------------------------------------------------- |
| 4609 | SUBROUTINE rc_gstom( icmp, compnam, lu, LAI_present, lai, glrad, sinphi, t, rh, diffc, & |
| 4610 | gstom, & |
| 4611 | p) |
| 4612 | |
| 4613 | ! input/output variables: |
| 4614 | INTEGER(iwp), INTENT(in) :: icmp ! component index |
| 4615 | CHARACTER(len=*), INTENT(in) :: compnam ! component name |
| 4616 | INTEGER(iwp), INTENT(in) :: lu ! land use type , lu = 1,...,nlu |
| 4617 | LOGICAL, INTENT(in) :: LAI_present |
| 4618 | REAL(kind=wp), INTENT(in) :: lai ! one-sided leaf area index |
| 4619 | REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W/m2) |
| 4620 | REAL(kind=wp), INTENT(in) :: sinphi ! sin of solar elevation angle |
| 4621 | REAL(kind=wp), INTENT(in) :: t ! temperature (C) |
| 4622 | REAL(kind=wp), INTENT(in) :: rh ! relative humidity (%) |
| 4623 | REAL(kind=wp), INTENT(in) :: diffc ! diffusion coefficient of the gas involved |
| 4624 | REAL(kind=wp), INTENT(out):: gstom ! stomatal conductance (m/s) |
| 4625 | REAL(kind=wp), OPTIONAL,INTENT(in) :: p ! pressure (Pa) |
| 4626 | |
| 4627 | |
| 4628 | ! Local variables |
| 4629 | REAL(kind=wp) :: vpd ! vapour pressure deficit (kPa) |
| 4630 | |
| 4631 | REAL(kind=wp), PARAMETER :: dO3 = 0.13e-4 ! diffusion coefficient of ozon (m2/s) |
| 4632 | |
| 4633 | |
| 4634 | SELECT CASE(trim(compnam)) |
| 4635 | |
| 4636 | !CASE('HNO3','N2O5','NO3','H2O2') this routine is not CALLed for HNO3 |
| 4637 | |
| 4638 | CASE('NO','CO') |
| 4639 | ! for NO stomatal uptake is neglected: |
| 4640 | gstom = 0.0 |
| 4641 | |
| 4642 | CASE('NO2','O3','SO2','NH3') |
| 4643 | |
| 4644 | ! IF vegetation present: |
| 4645 | IF (LAI_present) then |
| 4646 | |
| 4647 | IF (glrad .gt. 0.0) then |
| 4648 | CALL rc_get_vpd(t,rh,vpd) |
| 4649 | CALL rc_gstom_emb( lu, glrad, t, vpd, LAI_present, lai, sinphi, gstom, p ) |
| 4650 | gstom = gstom*diffc/dO3 ! Gstom of Emberson is derived for ozone |
| 4651 | ELSE |
| 4652 | gstom = 0.0 |
| 4653 | ENDIF |
| 4654 | ELSE |
| 4655 | ! no vegetation; zero conductance (infinite resistance): |
| 4656 | gstom = 0.0 |
| 4657 | ENDIF |
| 4658 | |
| 4659 | CASE default |
| 4660 | message_string = 'Component '// trim(compnam) // ' not supported' |
| 4661 | CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 ) |
| 4662 | END SELECT |
| 4663 | |
| 4664 | END SUBROUTINE rc_gstom |
| 4665 | |
| 4666 | |
| 4667 | |
| 4668 | !------------------------------------------------------------------- |
| 4669 | ! rc_gstom_emb: stomatal conductance according to Emberson |
| 4670 | !------------------------------------------------------------------- |
| 4671 | SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, LAI_present, lai, sinp, & |
| 4672 | Gsto, p ) |
| 4673 | ! History |
| 4674 | ! Original code from Lotos-Euros, TNO, M. Schaap |
| 4675 | ! 2009-08, M.C. van Zanten, Rivm |
| 4676 | ! Updated and extended. |
| 4677 | ! 2009-09, Arjo Segers, TNO |
| 4678 | ! Limitted temperature influence to range to avoid |
| 4679 | ! floating point exceptions. |
| 4680 | ! |
| 4681 | ! Method |
| 4682 | ! |
| 4683 | ! Code based on Emberson et al, 2000, Env. Poll., 403-413 |
| 4684 | ! Notation conform Unified EMEP Model Description Part 1, ch 8 |
| 4685 | ! |
| 4686 | ! In the calculation of F_light the modification of L. Zhang 2001, AE to the PARshade and PARsun |
| 4687 | ! parametrizations of Norman 1982 are applied |
| 4688 | ! |
| 4689 | ! F_phen and F_SWP are set to 1 |
| 4690 | ! |
| 4691 | ! Land use types DEPAC versus Emberson (Table 5.1, EMEP model description) |
| 4692 | ! DEPAC Emberson |
| 4693 | ! 1 = grass GR = grassland |
| 4694 | ! 2 = arable land TC = temperate crops ( LAI according to RC = rootcrops) |
| 4695 | ! 3 = permanent crops TC = temperate crops ( LAI according to RC = rootcrops) |
| 4696 | ! 4 = coniferous forest CF = tempareate/boREAL(kind=wp) coniferous forest |
| 4697 | ! 5 = deciduous forest DF = temperate/boREAL(kind=wp) deciduous forest |
| 4698 | ! 6 = water W = water |
| 4699 | ! 7 = urban U = urban |
| 4700 | ! 8 = other GR = grassland |
| 4701 | ! 9 = desert DE = desert |
| 4702 | ! |
| 4703 | |
| 4704 | ! Emberson specific declarations |
| 4705 | |
| 4706 | ! Input/output variables: |
| 4707 | INTEGER(iwp), INTENT(in) :: lu ! land use type, lu = 1,...,nlu |
| 4708 | REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W/m2) |
| 4709 | REAL(kind=wp), INTENT(in) :: T ! temperature (C) |
| 4710 | REAL(kind=wp), INTENT(in) :: vpd ! vapour pressure deficit (kPa) |
| 4711 | LOGICAL, INTENT(in) :: LAI_present |
| 4712 | REAL(kind=wp), INTENT(in) :: lai ! one-sided leaf area index |
| 4713 | REAL(kind=wp), INTENT(in) :: sinp ! sin of solar elevation angle |
| 4714 | REAL(kind=wp), INTENT(out):: Gsto ! stomatal conductance (m/s) |
| 4715 | REAL(kind=wp), OPTIONAL, INTENT(in) :: p ! pressure (Pa) |
| 4716 | |
| 4717 | ! Local variables: |
| 4718 | REAL(kind=wp) :: F_light, F_phen, F_temp, F_vpd, F_swp, bt, F_env |
| 4719 | REAL(kind=wp) :: pardir, pardiff, parshade, parsun, LAIsun, LAIshade, sinphi |
| 4720 | REAL(kind=wp) :: pres |
| 4721 | REAL(kind=wp), PARAMETER :: p_sealevel = 1.01325e5 ! Pa |
| 4722 | |
| 4723 | |
| 4724 | ! Check whether vegetation is present: |
| 4725 | IF (LAI_present) then |
| 4726 | |
| 4727 | ! calculation of correction factors for stomatal conductance |
| 4728 | IF (sinp .le. 0.0) then |
| 4729 | sinphi = 0.0001 |
| 4730 | ELSE |
| 4731 | sinphi = sinp |
| 4732 | END IF |
| 4733 | |
| 4734 | ! ratio between actual and sea-level pressure is used |
| 4735 | ! to correct for height in the computation of par ; |
| 4736 | ! should not exceed sea-level pressure therefore ... |
| 4737 | IF ( present(p) ) then |
| 4738 | pres = min( p, p_sealevel ) |
| 4739 | ELSE |
| 4740 | pres = p_sealevel |
| 4741 | ENDIF |
| 4742 | |
| 4743 | ! Direct and diffuse par, Photoactive (=visible) radiation: |
| 4744 | CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff ) |
| 4745 | |
| 4746 | ! par for shaded leaves (canopy averaged): |
| 4747 | parshade = pardiff*exp(-0.5*LAI**0.7)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Norman, 1982 |
| 4748 | IF (glrad .gt. 200 .and. LAI .gt. 2.5) then |
| 4749 | parshade = pardiff*exp(-0.5*LAI**0.8)+0.07*pardir*(1.1-0.1*LAI)*exp(-sinphi) ! Zhang et al., 2001 |
| 4750 | END IF |
| 4751 | |
| 4752 | ! par for sunlit leaves (canopy averaged): |
| 4753 | ! alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5 |
| 4754 | parsun = pardir*0.5/sinphi + parshade ! Norman, 1982 |
| 4755 | IF (glrad .gt. 200 .and. LAI .gt. 2.5) then |
| 4756 | parsun = pardir**0.8*0.5/sinphi + parshade ! Zhang et al., 2001 |
| 4757 | END IF |
| 4758 | |
| 4759 | ! leaf area index for sunlit and shaded leaves: |
| 4760 | IF (sinphi .gt. 0) then |
| 4761 | LAIsun = 2*sinphi*(1-exp(-0.5*LAI/sinphi )) |
| 4762 | LAIshade = LAI -LAIsun |
| 4763 | ELSE |
| 4764 | LAIsun = 0 |
| 4765 | LAIshade = LAI |
| 4766 | END IF |
| 4767 | |
| 4768 | F_light = (LAIsun*(1 - exp(-1.*alpha(lu)*parsun)) + LAIshade*(1 - exp(-1.*alpha(lu)*parshade)))/LAI |
| 4769 | |
| 4770 | F_light = max(F_light,F_min(lu)) |
| 4771 | |
| 4772 | ! temperature influence; only non-zero within range [Tmin,Tmax]: |
| 4773 | IF ( (Tmin(lu) < T) .and. (T < Tmax(lu)) ) then |
| 4774 | BT = (Tmax(lu)-Topt(lu))/(Topt(lu)-Tmin(lu)) |
| 4775 | F_temp = ((T-Tmin(lu))/(Topt(lu)-Tmin(lu))) * ((Tmax(lu)-T)/(Tmax(lu)-Topt(lu)))**BT |
| 4776 | ELSE |
| 4777 | F_temp = 0.0 |
| 4778 | END IF |
| 4779 | F_temp = max(F_temp,F_min(lu)) |
| 4780 | |
| 4781 | ! vapour pressure deficit influence |
| 4782 | F_vpd = min(1.,((1.-F_min(lu))*(vpd_min(lu)-vpd)/(vpd_min(lu)-vpd_max(lu)) + F_min(lu) )) |
| 4783 | F_vpd = max(F_vpd,F_min(lu)) |
| 4784 | |
| 4785 | F_swp = 1. |
| 4786 | |
| 4787 | ! influence of phenology on stom. conductance |
| 4788 | ! ignored for now in DEPAC since influence of F_phen on lu classes in use is negligible. |
| 4789 | ! When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore |
| 4790 | F_phen = 1. |
| 4791 | |
| 4792 | ! evaluate total stomatal conductance |
| 4793 | F_env = F_temp*F_vpd*F_swp |
| 4794 | F_env = max(F_env,F_min(lu)) |
| 4795 | gsto = G_max(lu) * F_light * F_phen * F_env |
| 4796 | |
| 4797 | ! gstom expressed per m2 leafarea; |
| 4798 | ! this is converted with LAI to m2 surface. |
| 4799 | Gsto = lai*gsto ! in m/s |
| 4800 | |
| 4801 | ELSE |
| 4802 | GSto = 0.0 |
| 4803 | ENDIF |
| 4804 | |
| 4805 | END SUBROUTINE rc_gstom_emb |
| 4806 | |
| 4807 | |
| 4808 | |
| 4809 | !------------------------------------------------------------------- |
| 4810 | ! par_dir_diff |
| 4811 | !------------------------------------------------------------------- |
| 4812 | SUBROUTINE par_dir_diff(glrad,sinphi,pres,pres_0,par_dir,par_diff) |
| 4813 | ! |
| 4814 | ! Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and |
| 4815 | ! diffuse, visible and near-infrared components. Agric. Forest Meteorol. |
| 4816 | ! 34, 205-213. |
| 4817 | ! |
| 4818 | ! From a SUBROUTINE obtained from Leiming Zhang (via Willem Asman), |
| 4819 | ! MeteoroLOGICAL Service of Canada |
| 4820 | ! e-mail: leiming.zhang@ec.gc.ca |
| 4821 | ! |
| 4822 | ! Leiming uses solar irradiance. This should be equal to global radiation and |
| 4823 | ! Willem Asman set it to global radiation |
| 4824 | |
| 4825 | REAL(kind=wp), INTENT(in) :: glrad ! global radiation (W m-2) |
| 4826 | REAL(kind=wp), INTENT(in) :: sinphi ! sine of the solar elevation |
| 4827 | REAL(kind=wp), INTENT(in) :: pres ! actual pressure (to correct for height) (Pa) |
| 4828 | REAL(kind=wp), INTENT(in) :: pres_0 ! pressure at sea level (Pa) |
| 4829 | REAL(kind=wp), INTENT(out) :: par_dir ! PAR direct : visible (photoactive) direct beam radiation (W m-2) |
| 4830 | REAL(kind=wp), INTENT(out) :: par_diff ! PAR diffuse: visible (photoactive) diffuse radiation (W m-2) |
| 4831 | |
| 4832 | ! fn = near-infrared direct beam fraction (DIMENSIONless) |
| 4833 | ! fv = PAR direct beam fraction (DIMENSIONless) |
| 4834 | ! ratio = ratio measured to potential solar radiation (DIMENSIONless) |
| 4835 | ! rdm = potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky |
| 4836 | ! rdn = potential diffuse near-infrared radiation (W m-2) |
| 4837 | ! rdu = visible (PAR) direct beam radiation (W m-2) |
| 4838 | ! rdv = potential visible (PAR) diffuse radiation (W m-2) |
| 4839 | ! rn = near-infrared radiation (W m-2) |
| 4840 | ! rv = visible radiation (W m-2) |
| 4841 | ! ww = water absorption in the near infrared for 10 mm of precipitable water |
| 4842 | |
| 4843 | REAL(kind=wp) :: rdu,rdv,ww,rdm,rdn,rv,rn,ratio,sv,fv |
| 4844 | |
| 4845 | ! Calculate visible (PAR) direct beam radiation |
| 4846 | ! 600 W m-2 represents average amount of PAR (400-700 nm wavelength) |
| 4847 | ! at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2) |
| 4848 | rdu=600.*exp(-0.185*(pres/pres_0)/sinphi)*sinphi |
| 4849 | |
| 4850 | ! Calculate potential visible diffuse radiation |
| 4851 | rdv=0.4*(600.- rdu)*sinphi |
| 4852 | |
| 4853 | ! Calculate the water absorption in the-near infrared |
| 4854 | ww=1320*10**( -1.195+0.4459*log10(1./sinphi)-0.0345*(log10(1./sinphi))**2 ) |
| 4855 | |
| 4856 | ! Calculate potential direct beam near-infrared radiation |
| 4857 | rdm=(720.*exp(-0.06*(pres/pres_0)/sinphi)-ww)*sinphi !720 = solar constant - 600 |
| 4858 | |
| 4859 | ! Calculate potential diffuse near-infrared radiation |
| 4860 | rdn=0.6*(720-rdm-ww)*sinphi |
| 4861 | |
| 4862 | ! Compute visible and near-infrared radiation |
| 4863 | rv=max(0.1,rdu+rdv) |
| 4864 | rn=max(0.01,rdm+rdn) |
| 4865 | |
| 4866 | ! Compute ratio between input global radiation and total radiation computed here |
| 4867 | ratio=min(0.9,glrad/(rv+rn)) |
| 4868 | |
| 4869 | ! Calculate total visible radiation |
| 4870 | sv=ratio*rv |
| 4871 | |
| 4872 | ! Calculate fraction of PAR in the direct beam |
| 4873 | fv=min(0.99, (0.9-ratio)/0.7) ! help variable |
| 4874 | fv=max(0.01,rdu/rv*(1.0-fv**0.6667)) ! fraction of PAR in the direct beam |
| 4875 | |
| 4876 | ! Compute direct and diffuse parts of PAR |
| 4877 | par_dir=fv*sv |
| 4878 | par_diff=sv-par_dir |
| 4879 | |
| 4880 | END SUBROUTINE par_dir_diff |
| 4881 | |
| 4882 | !------------------------------------------------------------------- |
| 4883 | ! rc_get_vpd: get vapour pressure deficit (kPa) |
| 4884 | !------------------------------------------------------------------- |
| 4885 | SUBROUTINE rc_get_vpd(temp, relh,vpd) |
| 4886 | |
| 4887 | IMPLICIT NONE |
| 4888 | |
| 4889 | ! Input/output variables: |
| 4890 | REAL(kind=wp) , INTENT(in) :: temp ! temperature (C) |
| 4891 | REAL(kind=wp) , INTENT(in) :: relh ! relative humidity (%) |
| 4892 | REAL(kind=wp) , INTENT(out) :: vpd ! vapour pressure deficit (kPa) |
| 4893 | |
| 4894 | ! Local variables: |
| 4895 | REAL(kind=wp) :: esat |
| 4896 | |
| 4897 | ! fit PARAMETERs: |
| 4898 | REAL(kind=wp), PARAMETER :: a1 = 6.113718e-1 |
| 4899 | REAL(kind=wp), PARAMETER :: a2 = 4.43839e-2 |
| 4900 | REAL(kind=wp), PARAMETER :: a3 = 1.39817e-3 |
| 4901 | REAL(kind=wp), PARAMETER :: a4 = 2.9295e-5 |
| 4902 | REAL(kind=wp), PARAMETER :: a5 = 2.16e-7 |
| 4903 | REAL(kind=wp), PARAMETER :: a6 = 3.0e-9 |
| 4904 | |
| 4905 | ! esat is saturation vapour pressure (kPa) at temp(C) following monteith(1973) |
| 4906 | esat = a1 + a2*temp + a3*temp**2 + a4*temp**3 + a5*temp**4 + a6*temp**5 |
| 4907 | vpd = esat * (1-relh/100) |
| 4908 | |
| 4909 | END SUBROUTINE rc_get_vpd |
| 4910 | |
| 4911 | |
| 4912 | |
| 4913 | !------------------------------------------------------------------- |
| 4914 | ! rc_gsoil_eff: compute effective soil conductance |
| 4915 | !------------------------------------------------------------------- |
| 4916 | SUBROUTINE rc_gsoil_eff(icmp,lu,sai,ust,nwet,t,gsoil_eff) |
| 4917 | |
| 4918 | ! Input/output variables: |
| 4919 | INTEGER(iwp) , INTENT(in) :: icmp ! component index |
| 4920 | INTEGER(iwp) , INTENT(in) :: lu ! land use type, lu = 1,...,nlu |
| 4921 | REAL(kind=wp) , INTENT(in) :: sai ! surface area index |
| 4922 | REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) |
| 4923 | INTEGER(iwp) , INTENT(in) :: nwet ! index for wetness |
| 4924 | ! nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow |
| 4925 | ! N.B. this routine cannot be CALLed with nwet = 9, |
| 4926 | ! nwet = 9 should be handled outside this routine. |
| 4927 | REAL(kind=wp) , INTENT(in) :: t ! temperature (C) |
| 4928 | REAL(kind=wp) , INTENT(out) :: gsoil_eff ! effective soil conductance (m/s) |
| 4929 | |
| 4930 | ! local variables: |
| 4931 | REAL(kind=wp) :: rinc ! in canopy resistance (s/m) |
| 4932 | REAL(kind=wp) :: rsoil_eff ! effective soil resistance (s/m) |
| 4933 | |
| 4934 | |
| 4935 | |
| 4936 | ! Soil resistance (numbers matched with lu_classes and component numbers) |
| 4937 | REAL(kind=wp), PARAMETER :: rsoil(nlu_dep,ncmp) = reshape( (/ & |
| 4938 | ! grs ara crp cnf dec wat urb oth des ice sav trf wai med sem |
| 4939 | 1000., 200., 200., 200., 200., 2000., 400., 1000., 2000., 2000., 1000., 200., 2000., 200., 400., & ! O3 |
| 4940 | 1000., 1000., 1000., 1000., 1000., 10., 1000., 1000., 1000., 500., 1000., 1000., 10., 1000., 1000., & ! SO2 |
| 4941 | 1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., & ! NO2 |
| 4942 | -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! NO |
| 4943 | 100., 100., 100., 100., 100., 10., 100., 100., 100., 1000., 100., 100., 10., 100., 100., & ! NH3 |
| 4944 | -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., & ! CO |
| 4945 | -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! NO3 |
| 4946 | -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! HNO3 |
| 4947 | -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., & ! N2O5 |
| 4948 | -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),& ! H2O2 |
| 4949 | (/nlu_dep,ncmp/) ) |
| 4950 | |
| 4951 | ! o3 so2 no2 no nh3 co no3 hno3 n2o5 h2o2 |
| 4952 | |
| 4953 | |
| 4954 | REAL(kind=wp), PARAMETER :: rsoil_wet(ncmp) = (/2000., 10., 2000., -999., 10., -999., -999., -999., -999., -999./) |
| 4955 | REAL(kind=wp), PARAMETER :: rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./) |
| 4956 | |
| 4957 | |
| 4958 | |
| 4959 | ! Compute in canopy (in crop) resistance: |
| 4960 | CALL rc_rinc(lu,sai,ust,rinc) |
| 4961 | |
| 4962 | ! Check for missing deposition path: |
| 4963 | IF (missing(rinc)) then |
| 4964 | rsoil_eff = -9999. |
| 4965 | ELSE |
| 4966 | ! Frozen soil (temperature below 0): |
| 4967 | IF (t .lt. 0.0) then |
| 4968 | IF (missing(rsoil_frozen(icmp))) then |
| 4969 | rsoil_eff = -9999. |
| 4970 | ELSE |
| 4971 | rsoil_eff = rsoil_frozen(icmp) + rinc |
| 4972 | ENDIF |
| 4973 | ELSE |
| 4974 | ! Non-frozen soil; dry: |
| 4975 | IF (nwet .eq. 0) then |
| 4976 | IF (missing(rsoil(lu,icmp))) then |
| 4977 | rsoil_eff = -9999. |
| 4978 | ELSE |
| 4979 | rsoil_eff = rsoil(lu,icmp) + rinc |
| 4980 | ENDIF |
| 4981 | |
| 4982 | ! Non-frozen soil; wet: |
| 4983 | ELSEIF (nwet .eq. 1) then |
| 4984 | IF (missing(rsoil_wet(icmp))) then |
| 4985 | rsoil_eff = -9999. |
| 4986 | ELSE |
| 4987 | rsoil_eff = rsoil_wet(icmp) + rinc |
| 4988 | ENDIF |
| 4989 | ELSE |
| 4990 | message_string = 'nwet can only be 0 or 1' |
| 4991 | CALL message( 'rc_gsoil_eff', 'CM0460', 1, 2, 0, 6, 0 ) |
| 4992 | ENDIF |
| 4993 | ENDIF |
| 4994 | ENDIF |
| 4995 | |
| 4996 | ! Compute conductance: |
| 4997 | IF (rsoil_eff .gt. 0.0) then |
| 4998 | gsoil_eff = 1./rsoil_eff |
| 4999 | ELSE |
| 5000 | gsoil_eff = 0.0 |
| 5001 | ENDIF |
| 5002 | |
| 5003 | END SUBROUTINE rc_gsoil_eff |
| 5004 | |
| 5005 | |
| 5006 | |
| 5007 | !------------------------------------------------------------------- |
| 5008 | ! rc_rinc: compute in canopy (or in crop) resistance |
| 5009 | ! van Pul and Jacobs, 1993, BLM |
| 5010 | !------------------------------------------------------------------- |
| 5011 | SUBROUTINE rc_rinc(lu,sai,ust,rinc) |
| 5012 | |
| 5013 | |
| 5014 | ! Input/output variables: |
| 5015 | INTEGER(iwp) , INTENT(in) :: lu ! land use class, lu = 1, ..., nlu |
| 5016 | REAL(kind=wp) , INTENT(in) :: sai ! surface area index |
| 5017 | REAL(kind=wp) , INTENT(in) :: ust ! friction velocity (m/s) |
| 5018 | REAL(kind=wp) , INTENT(out) :: rinc ! in canopy resistance (s/m) |
| 5019 | |
| 5020 | |
| 5021 | ! b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 IF not applicable) |
| 5022 | ! h = vegetation height (m) grass arabl crops conIF decid water urba othr desr ice sav trf wai med semi |
| 5023 | REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999, 14, 14 /) |
| 5024 | REAL(kind=wp), DIMENSION(nlu_dep), PARAMETER :: h = (/ -999, 1, 1, 20, 20, -999, -999, -999, -999, -999, -999, 20, -999, 1 , 1 /) |
| 5025 | |
| 5026 | ! Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0: |
| 5027 | IF (b(lu) .gt. 0.0) then |
| 5028 | |
| 5029 | ! Check for u* > 0 (otherwise denominator = 0): |
| 5030 | IF (ust .gt. 0.0) then |
| 5031 | rinc = b(lu)*h(lu)*sai/ust |
| 5032 | ELSE |
| 5033 | rinc = 1000.0 |
| 5034 | ENDIF |
| 5035 | ELSE |
| 5036 | IF (lu .eq. ilu_grass .or. lu .eq. ilu_other ) then |
| 5037 | rinc = -999. ! no deposition path for grass, other, and semi-natural |
| 5038 | ELSE |
| 5039 | rinc = 0.0 ! no in-canopy resistance |
| 5040 | ENDIF |
| 5041 | ENDIF |
| 5042 | |
| 5043 | END SUBROUTINE rc_rinc |
| 5044 | |
| 5045 | |
| 5046 | |
| 5047 | !------------------------------------------------------------------- |
| 5048 | ! rc_rctot: compute total canopy (or surface) resistance Rc |
| 5049 | !------------------------------------------------------------------- |
| 5050 | SUBROUTINE rc_rctot(gstom,gsoil_eff,gw,gc_tot,rc_tot) |
| 5051 | |
| 5052 | ! Input/output variables: |
| 5053 | REAL(kind=wp), INTENT(in) :: gstom ! stomatal conductance (s/m) |
| 5054 | REAL(kind=wp), INTENT(in) :: gsoil_eff ! effective soil conductance (s/m) |
| 5055 | REAL(kind=wp), INTENT(in) :: gw ! external leaf conductance (s/m) |
| 5056 | REAL(kind=wp), INTENT(out) :: gc_tot ! total canopy conductance (m/s) |
| 5057 | REAL(kind=wp), INTENT(out) :: rc_tot ! total canopy resistance Rc (s/m) |
| 5058 | |
| 5059 | ! Total conductance: |
| 5060 | gc_tot = gstom + gsoil_eff + gw |
| 5061 | |
| 5062 | ! Total resistance (note: gw can be negative, but no total emission allowed here): |
| 5063 | IF (gc_tot .le. 0.0 .or. gw .lt. 0.0) then |
| 5064 | rc_tot = -9999. |
| 5065 | ELSE |
| 5066 | rc_tot = 1./gc_tot |
| 5067 | ENDIF |
| 5068 | |
| 5069 | END SUBROUTINE rc_rctot |
| 5070 | |
| 5071 | |
| 5072 | |
| 5073 | !------------------------------------------------------------------- |
| 5074 | ! rc_comp_point_rc_eff: calculate the effective resistance Rc |
| 5075 | ! based on one or more compensation points |
| 5076 | !------------------------------------------------------------------- |
| 5077 | ! |
| 5078 | ! NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173. |
| 5079 | ! Sutton 1998 AE 473-480) |
| 5080 | ! |
| 5081 | ! Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine. |
| 5082 | ! FS 2009-01-29: variable names made consistent with DEPAC |
| 5083 | ! FS 2009-03-04: use total compensation point |
| 5084 | ! |
| 5085 | ! C: with total compensation point ! D: approximation of C |
| 5086 | ! ! with classical approach |
| 5087 | ! zr --------- Catm ! zr --------- Catm |
| 5088 | ! | ! | |
| 5089 | ! Ra ! Ra |
| 5090 | ! | ! | |
| 5091 | ! Rb ! Rb |
| 5092 | ! | ! | |
| 5093 | ! z0 --------- Cc ! z0 --------- Cc |
| 5094 | ! | ! | |
| 5095 | ! Rc ! Rc_eff |
| 5096 | ! | ! | |
| 5097 | ! --------- Ccomp_tot ! --------- C=0 |
| 5098 | ! |
| 5099 | ! |
| 5100 | ! The effective Rc is defined such that instead of using |
| 5101 | ! |
| 5102 | ! F = -vd*[Catm - Ccomp_tot] (1) |
| 5103 | ! |
| 5104 | ! we can use the 'normal' flux formula |
| 5105 | ! |
| 5106 | ! F = -vd'*Catm, (2) |
| 5107 | ! |
| 5108 | ! with vd' = 1/(Ra + Rb + Rc') (3) |
| 5109 | ! |
| 5110 | ! and Rc' the effective Rc (rc_eff). |
| 5111 | ! (Catm - Ccomp_tot) |
| 5112 | ! vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------ |
| 5113 | ! Catm |
| 5114 | ! |
| 5115 | ! (Catm - Ccomp_tot) |
| 5116 | ! 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------ |
| 5117 | ! Catm |
| 5118 | ! |
| 5119 | ! Catm |
| 5120 | ! (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------ |
| 5121 | ! (Catm - Ccomp_tot) |
| 5122 | ! |
| 5123 | ! Catm |
| 5124 | ! Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb |
| 5125 | ! (Catm - Ccomp_tot) |
| 5126 | ! |
| 5127 | ! Catm Catm |
| 5128 | ! Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------ |
| 5129 | ! (Catm - Ccomp_tot) (Catm - Ccomp_tot) |
| 5130 | ! |
| 5131 | ! Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot) |
| 5132 | ! |
| 5133 | ! This is not the most efficient way to DO this; |
| 5134 | ! in the current LE version, (1) is used directly in the CALLing routine |
| 5135 | ! and this routine is not used anymore. |
| 5136 | ! |
| 5137 | ! ------------------------------------------------------------------------------------------- |
| 5138 | ! In fact it is the question IF this correct; note the difference in differential equations |
| 5139 | ! |
| 5140 | ! dCatm ! dCatm |
| 5141 | ! H ----- = F = -vd'*Catm, ! H ----- = F = -vd*[Catm - Ccomp_tot], |
| 5142 | ! dt ! dt |
| 5143 | ! |
| 5144 | ! where in the left colum vd' is a function of Catm and not a constant anymore. |
| 5145 | ! |
| 5146 | ! Another problem is that IF Catm -> 0, we end up with an infinite value of the exchange |
| 5147 | ! velocity vd. |
| 5148 | ! ------------------------------------------------------------------------------------------- |
| 5149 | |
| 5150 | SUBROUTINE rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) |
| 5151 | |
| 5152 | IMPLICIT NONE |
| 5153 | |
| 5154 | ! Input/output variables: |
| 5155 | REAL(kind=wp), INTENT(in) :: ccomp_tot ! total compensation point (weighed average of separate compensation points) (ug/m3) |
| 5156 | REAL(kind=wp), INTENT(in) :: catm ! atmospheric concentration (ug/m3) |
| 5157 | REAL(kind=wp), INTENT(in) :: ra ! aerodynamic resistance (s/m) |
| 5158 | REAL(kind=wp), INTENT(in) :: rb ! boundary layer resistance (s/m) |
| 5159 | REAL(kind=wp), INTENT(in) :: rc_tot ! total canopy resistance (s/m) |
| 5160 | REAL(kind=wp), INTENT(out) :: rc_eff ! effective total canopy resistance (s/m) |
| 5161 | |
| 5162 | ! Compute effective resistance: |
| 5163 | IF ( ccomp_tot == 0.0 ) then |
| 5164 | ! trace with no compensiation point ( or compensation point equal to zero) |
| 5165 | rc_eff = rc_tot |
| 5166 | |
| 5167 | ELSE IF ( ccomp_tot > 0 .and. ( abs( catm - ccomp_tot ) < 1.e-8 ) ) then |
| 5168 | ! surface concentration (almoast) equal to atmospheric concentration |
| 5169 | ! no exchange between surface and atmosphere, infinite RC --> vd=0 |
| 5170 | rc_eff = 9999999999. |
| 5171 | |
| 5172 | ELSE IF ( ccomp_tot > 0 ) then |
| 5173 | ! compensation point available, calculate effective restisance |
| 5174 | rc_eff = ((ra + rb)*ccomp_tot + rc_tot*catm)/(catm-ccomp_tot) |
| 5175 | |
| 5176 | ELSE |
| 5177 | rc_eff = -999. |
| 5178 | message_string = 'This should not be possible, check ccomp_tot' |
| 5179 | CALL message( 'rc_comp_point_rc_eff', 'CM0461', 1, 2, 0, 6, 0 ) |
| 5180 | ENDIF |
| 5181 | |
| 5182 | return |
| 5183 | END SUBROUTINE rc_comp_point_rc_eff |
| 5184 | |
| 5185 | |
| 5186 | |
| 5187 | !------------------------------------------------------------------- |
| 5188 | ! missing: check for data that correspond with a missing deposition path |
| 5189 | ! this data is represented by -999 |
| 5190 | !------------------------------------------------------------------- |
| 5191 | |
| 5192 | LOGICAL function missing(x) |
| 5193 | |
| 5194 | REAL(kind=wp), INTENT(in) :: x |
| 5195 | |
| 5196 | ! bandwidth for checking (in)equalities of floats |
| 5197 | REAL(kind=wp), PARAMETER :: EPS = 1.0e-5 |
| 5198 | |
| 5199 | missing = (abs(x + 999.) .le. EPS) |
| 5200 | |
| 5201 | END function missing |
| 5202 | |
| 5203 | |
| 5204 | |
| 5205 | ELEMENTAL FUNCTION sedimentation_velocity(rhopart, partsize, slipcor, visc) RESULT (vs) |
| 5206 | |
| 5207 | |
| 5208 | IMPLICIT NONE |
| 5209 | |
| 5210 | ! --- in/out --------------------------------- |
| 5211 | |
| 5212 | REAL(kind=wp), INTENT(in) :: rhopart ! kg/m3 |
| 5213 | REAL(kind=wp), INTENT(in) :: partsize ! m |
| 5214 | REAL(kind=wp), INTENT(in) :: slipcor ! m |
| 5215 | REAL(kind=wp), INTENT(in) :: visc ! viscosity |
| 5216 | REAL(kind=wp) :: vs |
| 5217 | |
| 5218 | ! acceleration of gravity: |
| 5219 | REAL(kind=wp), PARAMETER :: grav = 9.80665 ! m/s2 |
| 5220 | |
| 5221 | ! --- begin ---------------------------------- |
| 5222 | |
| 5223 | !viscosity & sedimentation velocity |
| 5224 | vs = rhopart * (partsize**2) * grav * slipcor / (18*visc) |
| 5225 | |
| 5226 | END FUNCTION sedimentation_velocity |
| 5227 | |
| 5228 | !------------------------------------------------------------------------ |
| 5229 | ! Boundary-layer deposition resistance following Zhang (2001) |
| 5230 | !------------------------------------------------------------------------ |
| 5231 | SUBROUTINE drydepo_aero_zhang_vd( vd, Rs, vs1, partsize, slipcor, & |
| 5232 | nwet, tsurf, dens1, viscos1, & |
| 5233 | luc, ftop_lu, ustar_lu) |
| 5234 | |
| 5235 | |
| 5236 | IMPLICIT NONE |
| 5237 | |
| 5238 | ! --- in/out --------------------------------- |
| 5239 | |
| 5240 | REAL(kind=wp), INTENT(out) :: vd ! deposition velocity (m/s) |
| 5241 | REAL(kind=wp), INTENT(out) :: Rs ! Sedimentaion resistance (s/m) |
| 5242 | REAL(kind=wp), INTENT(in) :: vs1 ! sedimentation velocity in lowest layer |
| 5243 | REAL(kind=wp), INTENT(in) :: partsize ! particle diameter (m) |
| 5244 | REAL(kind=wp), INTENT(in) :: slipcor ! slip correction factor |
| 5245 | REAL(kind=wp), INTENT(in) :: tsurf ! surface temperature (K) |
| 5246 | REAL(kind=wp), INTENT(in) :: dens1 ! air density (kg/m3) in lowest layer |
| 5247 | REAL(kind=wp), INTENT(in) :: viscos1 ! air viscosity in lowest layer |
| 5248 | REAL(kind=wp), INTENT(in) :: ftop_lu ! atmospheric resistnace Ra |
| 5249 | REAL(kind=wp), INTENT(in) :: ustar_lu ! friction velocity u* |
| 5250 | |
| 5251 | INTEGER(iwp), INTENT(in) :: nwet ! 1=rain, 9=snowcover |
| 5252 | INTEGER(iwp), INTENT(in) :: luc ! DEPAC LU |
| 5253 | |
| 5254 | |
| 5255 | ! --- const ---------------------------------- |
| 5256 | |
| 5257 | |
| 5258 | |
| 5259 | ! acceleration of gravity: |
| 5260 | REAL(kind=wp), PARAMETER :: grav = 9.80665 ! m/s2 |
| 5261 | |
| 5262 | REAL(kind=wp), PARAMETER :: beta = 2.0 |
| 5263 | REAL(kind=wp), PARAMETER :: epsilon0 = 3.0 |
| 5264 | REAL(kind=wp), PARAMETER :: kb = 1.38066e-23 |
| 5265 | REAL(kind=wp), PARAMETER :: pi = 3.141592654_wp !< PI |
| 5266 | |
| 5267 | REAL, PARAMETER :: alfa_lu(nlu_dep) = & |
| 5268 | (/1.2, 1.2, 1.2, 1.0, 1.0, 100.0, 1.5, 1.2, 50.0, 100.0, 1.2, 1.0, 100.0, 1.2, 50.0/) |
| 5269 | ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem |
| 5270 | REAL, PARAMETER :: gamma_lu(nlu_dep) = & |
| 5271 | (/0.54, 0.54, 0.54, 0.56, 0.56, 0.50, 0.56, 0.54, 0.58, 0.50, 0.54, 0.56, 0.50, 0.54, 0.54/) |
| 5272 | ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem |
| 5273 | REAL, PARAMETER ::A_lu(nlu_dep) = & |
| 5274 | (/3.0, 3.0, 2.0, 2.0, 7.0, -99., 10.0, 3.0, -99., -99., 3.0, 7.0, -99., 2.0, -99./) |
| 5275 | ! grass arabl crops conif decid water urba othr desr ice sav trf wai med sem |
| 5276 | |
| 5277 | |
| 5278 | ! --- local ---------------------------------- |
| 5279 | |
| 5280 | REAL(kind=wp) :: kinvisc, diff_part |
| 5281 | REAL(kind=wp) :: schmidt,stokes, Ebrown, Eimpac, Einterc, Reffic |
| 5282 | |
| 5283 | ! --- begin ---------------------------------- |
| 5284 | |
| 5285 | ! kinetic viscosity & diffusivity |
| 5286 | kinvisc = viscos1 / dens1 ! only needed at surface |
| 5287 | |
| 5288 | diff_part = kb * tsurf * slipcor / (3*pi*viscos1*partsize) |
| 5289 | |
| 5290 | ! Schmidt number |
| 5291 | schmidt = kinvisc / diff_part |
| 5292 | |
| 5293 | ! calculate collection efficiencie E |
| 5294 | Ebrown = Schmidt**(-gamma_lu(luc)) ! Brownian diffusion |
| 5295 | ! determine Stokes number, interception efficiency |
| 5296 | ! and sticking efficiency R (1 = no rebound) |
| 5297 | IF ( luc == ilu_ice .or. nwet.eq.9 .or. luc == ilu_water_sea .or. luc == ilu_water_inland ) THEN |
| 5298 | stokes=vs1*ustar_lu**2/(grav*kinvisc) |
| 5299 | Einterc=0.0 |
| 5300 | Reffic=1.0 |
| 5301 | ELSE IF ( luc == ilu_other .or. luc == ilu_desert ) THEN !tundra of desert |
| 5302 | stokes=vs1*ustar_lu**2/(grav*kinvisc) |
| 5303 | Einterc=0.0 |
| 5304 | Reffic=exp(-Stokes**0.5) |
| 5305 | ELSE |
| 5306 | stokes=vs1*ustar_lu/(grav*A_lu(luc)*1.e-3) |
| 5307 | Einterc=0.5*(partsize/(A_lu(luc)*1e-3))**2 |
| 5308 | Reffic=exp(-Stokes**0.5) |
| 5309 | END IF |
| 5310 | ! when surface is wet all particles DO not rebound: |
| 5311 | IF(nwet.eq.1) Reffic = 1.0 |
| 5312 | ! determine impaction efficiency: |
| 5313 | Eimpac = ( stokes / (alfa_lu(luc)+stokes) )**beta |
| 5314 | |
| 5315 | ! sedimentation resistance: |
| 5316 | Rs = 1.0 / ( epsilon0 * ustar_lu * (Ebrown+Eimpac+Einterc) * Reffic ) |
| 5317 | |
| 5318 | ! deposition velocity according to Seinfeld and Pandis (= SP, 2006; eq 19.7): |
| 5319 | ! |
| 5320 | ! 1 |
| 5321 | ! vd = ------------------ + vs |
| 5322 | ! Ra + Rs + Ra*Rs*vs |
| 5323 | ! |
| 5324 | ! where: Rs = Rb (in SP, 2006) |
| 5325 | ! |
| 5326 | ! |
| 5327 | vd = 1.0 / ( ftop_lu + Rs + ftop_lu*Rs*vs1) + vs1 |
| 5328 | |
| 5329 | |
| 5330 | |
| 5331 | END SUBROUTINE drydepo_aero_zhang_vd |
| 5332 | |
| 5333 | |
| 5334 | |
| 5335 | SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb ) |
| 5336 | |
| 5337 | ! Compute quasi-laminar boundary layer resistance as a function of landuse and tracer. |
| 5338 | ! |
| 5339 | ! Original EMEP formulation by (Simpson et al, 2003) is used. |
| 5340 | ! |
| 5341 | |
| 5342 | |
| 5343 | IMPLICIT NONE |
| 5344 | |
| 5345 | ! --- in/out --------------------------------- |
| 5346 | |
| 5347 | LOGICAL , INTENT(in) :: is_water |
| 5348 | REAL(kind=wp), INTENT(in) :: z0h ! roughness length for heat |
| 5349 | REAL(kind=wp), INTENT(in) :: ustar ! friction velocity |
| 5350 | REAL(kind=wp), INTENT(in) :: diffc ! coefficient of diffusivity |
| 5351 | REAL(kind=wp), INTENT(out) :: Rb ! boundary layer resistance |
| 5352 | |
| 5353 | ! --- const ---------------------------------- |
| 5354 | |
| 5355 | ! thermal diffusivity of dry air 20 C |
| 5356 | REAL(kind=wp), PARAMETER :: thk = 0.19e-4 ! http://en.wikipedia.org/wiki/Thermal_diffusivity (@ T=300K) |
| 5357 | REAL(kind=wp), PARAMETER :: kappa_stab = 0.35 ! von Karman constant |
| 5358 | |
| 5359 | ! --- begin --------------------------------- |
| 5360 | |
| 5361 | |
| 5362 | ! Use Simpson et al. (2003) |
| 5363 | !IF ( is_water ) THEN |
| 5364 | !Rb = 1.0_wp / (kappa_stab*max(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.01_wp,ustar)) |
| 5365 | !TODO: Check Rb over water calculation!!! |
| 5366 | ! Rb = 1.0_wp / (kappa_stab*max(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*max(0.1_wp,ustar)) |
| 5367 | !ELSE |
| 5368 | Rb = 5.0_wp / max(0.01_wp,ustar) * (thk/diffc)**0.67_wp |
| 5369 | !END IF |
| 5370 | |
| 5371 | END SUBROUTINE get_rb_cell |
| 5372 | |
| 5373 | |
| 5374 | !----------------------------------------------------------------- |
| 5375 | ! Compute water vapor partial pressure (e_w) |
| 5376 | ! given specific humidity Q [(kg water)/(kg air)]. |
| 5377 | ! |
| 5378 | ! Use that gas law for volume V with temperature T |
| 5379 | ! holds for the total mixture as well as the water part: |
| 5380 | ! |
| 5381 | ! R T / V = p_air / n_air = p_water / n_water |
| 5382 | ! |
| 5383 | ! thus: |
| 5384 | ! |
| 5385 | ! p_water = p_air n_water / n_air |
| 5386 | ! |
| 5387 | ! Use: |
| 5388 | ! n_air = m_air / xm_air |
| 5389 | ! [kg air] / [(kg air)/(mole air)] |
| 5390 | ! and: |
| 5391 | ! n_water = m_air * Q / xm_water |
| 5392 | ! [kg water] / [(kg water)/(mole water)] |
| 5393 | ! thus: |
| 5394 | ! p_water = p_air Q / (xm_water/xm_air) |
| 5395 | ! |
| 5396 | |
| 5397 | ELEMENTAL FUNCTION watervaporpartialpressure( Q, p ) result ( p_w ) |
| 5398 | |
| 5399 | |
| 5400 | ! --- in/out --------------------------------- |
| 5401 | |
| 5402 | REAL(kind=wp), INTENT(in) :: Q ! specific humidity [(kg water)/(kg air)] |
| 5403 | REAL(kind=wp), INTENT(in) :: p ! air pressure [Pa] |
| 5404 | REAL(kind=wp) :: p_w ! water vapor partial pressure [Pa] |
| 5405 | |
| 5406 | ! --- const ---------------------------------- |
| 5407 | |
| 5408 | ! mole mass ratio: |
| 5409 | REAL(kind=wp), PARAMETER :: eps = xm_H2O / xm_air ! ~ 0.622 |
| 5410 | |
| 5411 | ! --- begin ---------------------------------- |
| 5412 | |
| 5413 | ! partial pressure of water vapor: |
| 5414 | p_w = p * Q / eps |
| 5415 | |
| 5416 | END function watervaporpartialpressure |
| 5417 | |
| 5418 | |
| 5419 | |
| 5420 | !------------------------------------------------------------------ |
| 5421 | ! Saturation vapor pressure. |
| 5422 | ! From (Stull 1988, eq. 7.5.2d): |
| 5423 | ! |
| 5424 | ! e_sat = p0 exp( 17.67 * (T-273.16) / (T-29.66) ) [Pa] |
| 5425 | ! |
| 5426 | ! where: |
| 5427 | ! p0 = 611.2 [Pa] : reference pressure |
| 5428 | ! |
| 5429 | ! Arguments: |
| 5430 | ! T [K] : air temperature |
| 5431 | ! Result: |
| 5432 | ! e_sat_w [Pa] : saturation vapor pressure |
| 5433 | ! |
| 5434 | ! References: |
| 5435 | ! Roland B. Stull, 1988 |
| 5436 | ! An introduction to boundary layer meteorology. |
| 5437 | ! |
| 5438 | |
| 5439 | ELEMENTAL FUNCTION saturationvaporpressure( T ) result( e_sat_w ) |
| 5440 | |
| 5441 | |
| 5442 | ! --- in/out --------------------------------- |
| 5443 | |
| 5444 | REAL(kind=wp), INTENT(in) :: T ! temperature [K] |
| 5445 | REAL(kind=wp) :: e_sat_w ! saturation vapor pressure [Pa] |
| 5446 | |
| 5447 | ! --- const ---------------------------------- |
| 5448 | |
| 5449 | ! base pressure: |
| 5450 | REAL(kind=wp), PARAMETER :: p0 = 611.2 ! [Pa] |
| 5451 | |
| 5452 | ! --- begin ---------------------------------- |
| 5453 | |
| 5454 | ! saturation vapor pressure: |
| 5455 | e_sat_w = p0 * exp( 17.67 * (T-273.16) / (T-29.66) ) ! [Pa] |
| 5456 | |
| 5457 | END FUNCTION saturationvaporpressure |
| 5458 | |
| 5459 | |
| 5460 | |
| 5461 | !------------------------------------------------------------------------ |
| 5462 | ! Relative humidity RH [%] is by definition: |
| 5463 | ! |
| 5464 | ! e_w # water vapor partial pressure |
| 5465 | ! Rh = -------- * 100 |
| 5466 | ! e_sat_w # saturation vapor pressure |
| 5467 | ! |
| 5468 | ! Compute here from: |
| 5469 | ! Q specific humidity [(kg water)/(kg air)] |
| 5470 | ! T temperature [K] |
| 5471 | ! P pressure [Pa] |
| 5472 | ! |
| 5473 | |
| 5474 | ELEMENTAL FUNCTION relativehumidity_from_specifichumidity( Q, T, p ) result( Rh ) |
| 5475 | |
| 5476 | |
| 5477 | ! --- in/out --------------------------------- |
| 5478 | |
| 5479 | REAL(kind=wp), INTENT(in) :: Q ! specific humidity [(kg water)/(kg air)] |
| 5480 | REAL(kind=wp), INTENT(in) :: T ! temperature [K] |
| 5481 | REAL(kind=wp), INTENT(in) :: p ! air pressure [Pa] |
| 5482 | REAL(kind=wp) :: Rh ! relative humidity [%] |
| 5483 | |
| 5484 | ! --- begin ---------------------------------- |
| 5485 | |
| 5486 | ! relative humidity: |
| 5487 | Rh = watervaporpartialpressure( Q, p ) / saturationvaporpressure( T ) * 100.0 |
| 5488 | |
| 5489 | END FUNCTION relativehumidity_from_specifichumidity |
| 5490 | |
| 5491 | |
| 5492 | |